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Abstract 


A  frozen  orbit  is  an  orbit  whose  time  rate  of  change  of  the  argument  of  the 
periapsis  ) ,  the  eccentricity  (e),  the  semi  major  axis  (a),  or  the  angle  of  inclination  (0  is 
approximately  equal  to  zero.  Martian  frozen  orbits  are  known  to  exist  for  polar 
trajectories  with  altitudes  from  300  km  to  1000  km.  The  objective  of  this  study  was  to 
determine  if  other  regions  with  characteristics  similar  to  the  known  frozen  orbits  exist, 
taking  into  account  the  perturbative  effects  due  to  a  6  X  6  gravity  field  and  atmospheric 
drag. 

First,  the  geopotential  equation  was  derived  for  both  spherical  coordinates  and  the 
classical  orbital  elements.  Next,  a  model  for  the  atmospheric  drag  was  developed.  Using 
these  two  models,  a  Fortran  computer  model  named  ASAP  (Artificial  Satellite  Analysis 
Program)  was  analyzed  for  accuracy.  This  program  proved  to  be  highly  reliable,  and 
was  used  to  carry  out  further  analysis. 

Two  of  the  three  trajectories  planned  for  the  future  Mars  Geoscience/Climatology 
Orbiter  (MGCO)  are  frozen  orbits.  In  order  to  determine  the  characteristics  of  u>,  e,  a, 
and  i  of  a  frozen  orbit,  one  of  the  MGCO  frozen  orbits  was  examined  in  both  a  6  X  0 
and  a  6  X  6  gravity  field.  The  analysis  showed  that  the  above  orbital  elements  are  not 
periodic  over  one  orbital  period  (when  in  the  presence  of  a  6  X  6  gravity  field),  but 
they  are  bounded  over  one  axial  period.  ■  ^  » 

Since  the  greatest  change  over  an  orbital  period  is  in  the  argument  of  the  periapsis 
and  the  eccentricity  the  effect  of  driving  the  change  in  these  two  parameters  to 
approximately  zero  over  one  orbital  period  was  investigated.  Driving  the  change  in  u..  to 
zero  does  not  provide  the  desired  level  of  control  on  the  argument  of  the  periapsis. 
Driving  the  change  in  «  to  zero  can  only  be  accomplished  at  the  cost  of  relatively  high 


rates  of  change  in  u>  over  one  orbital  period.  An  orbit  was  found  in  which  the  change 
in  o>  and  in  e  over  one  orbital  period  were  both  equal  to  approximately  zero.  Again  the 
argument  of  the  periapsis  is  not  bounded,  but  rather  periodic. 


A  search  for  a  combination  of  orbital  elements  which  would  yield  a  zero  change 


over  one  orbital  period  for  all  four  of  the  above  orbital  elements  was  conducted  for  an 
eccentricity  of  0.3.  The  results  showed  no  such  orbit  exist;  regions  were  found  in 
which  the  change  in  3  out  of  the  4  orbital  elements  were  driven  to  zero  nr 
approximately  zero. 

Finally,  the  predominant  characteristics  of  the  elements  for  the  MGCO  frozen  orbit 
are  identified,  and  a  region  with  these  same  characteristics  was  found. 
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FROZEN  ORBIT  ANALYSIS  IN  THE  MARTIAN  SYSTEM 


I.  Introduction 


Background 

Mars  is  the  closest  planet  to  Earth  that  is  potentially  habitable  by  man;  however. 
Mars  is,  at  its  closest  approach  to  Earth,  approximately  78  million  kilometers  away. 
Even  though  the  trip  to  Mars  is  made  along  keplerian  trajectories  which  take  advantage 
of  the  Sun’s  gravity,  fuel  is  still  consumed  in  trajectory  correction  maneuvers.  When  a 
probe  arrives  at  Mars,  fuel  will  again  be  required  to  establish,  and  maintain  an  orbit 
about  the  planet.  The  size  of  the  probe  that  can  be  sent  to  Mars  is  dependent  on  the 
size  of  the  booster  used  to  get  the  probe  out  of  the  Earth’s  gravity  field.  Mission 
planners  must  make  use  of  boosters  currently  available  because  both  budget  and  time 
constraints  do  not  allow  for  a  booster  to  be  designed  for  a  specific  mission.  Therefore, 
the  size  of  the  payload  is  itself  a  constraint,  part  of  which  is  taken  up  in  mission 
required  fuel.  If  the  mission  profile  is  such  that  the  fuel  required  is  minimized,  then 
the  mission  duration  can  be  increased.  Having  the  capability  to  maintain  probes  in  orbit 
about  Mars  for  long  periods  will  increase  our  knowledge  of  Mars’  surface,  climatology, 
gravity  field,  magnetic  field,  and  the  interaction  of  the  magnetic  field  with  the  solar 
wind.  Such  a  probe  could  also  be  used  to  better  determine  what  and  where  Mars’ 
resources  are,  a  factor  that  may  be  critical  to  future  manned  missions  to  the  planet. 

One  method  of  minimizing  fuel  is  to  select  an  orbit  that  takes  advantage  of  Mars’ 
gravity  field  in  such  a  way  as  to  minimize  the  effects  of  atmospheric  drag  upon  the 
probe.  In  the  1960’s  H.  W.  West,  R.  T.  Clapp,  and  H.  Small  were  able  to  show  that  for 
Earth  there  exist  a  class  of  polar  orbits  with  non  zero  eccentricity,  and  whose  argument 
of  the  periapsis  is  over  the  south  pole,  such  that  the  line  of  apsides  does  not  rotate,  but 
rather  oscillates  about  its  initial  position.  These  orbits  were  called  "frozen"  because  of 


the  off  setting  effects  of  the  odd  and  even  zonal  harmonics  on  the  eccentricity  and  the 
argument  of  the  periapsis  yielding  orbits  whose  shape,  and  whose  orientation  of  the  line 
of  apsides  is  nearly  constant  over  time  (17:2).  Since  most  planets  are  oblate,  and  since 
atmospheric  drag  is  a  function  of  altitude  above  the  planet,  the  amount  of  atmospheric 
drag  experienced  will  be  less  over  the  poles.  This  implies  that  a  probe  in  an  orbit  that 
maintains  its  periapsis  over  a  polar  region  will  experience  less  drag,  and  hence  require 
less  fuel  consumption  to  remain  in  orbit.  For  Mars,  frozen  orbits  are  known  to  exist  for 
polar,  or  near  polar  orbits  with  altitudes  from  300  to  1000  km  (17:2). 


Definition  of  a  Frozen  Orbit 

This  thesis  defines  a  frozen  orbit  as  any  orbit  whose  time  rate  of  change  of  the 
argument  of  the  periapsis  (^),  the  eccentricity  (e),  the  semi  major  axis  (a),  or  the  angle 
of  inclination  (<)  is  equal  to  approximately  zero. 

Objective 

Given  the  perturbing  effects  of  the  zonal  and  sectoral  harmonics  up  to  and 
including  an  order  of  six,  and  the  perturbing  effects  of  atmospheric  drag,  this  thesis 
seeks  to  determine  other  regions  where  orbital  stabilities  similar  to  the  polar  frozen 
orbits  may  exist. 

Methodology 

First,  in  order  to  understand  the  relationship  that  exist  between  the  orbital 
elements  for  a  frozen  orbit,  a  known  Martian  frozen  orbit  will  be  examined.  From  the 
understanding  of  the  sensitivities  of  this  orbit  to  changes  in  the  orbital  elements, 
manipulations  of  the  orbital  elements  will  be  made  in  an  effort  to  find  other  stable 
regions.  This  thesis  will  only  consider  the  perturbing  effects  due  to  the  geopotential  and 
atmospheric  drag  upon  orbits  with  altitudes  from  approximately  200  km  to  20,000  km 
(Martian  geosynchronous).  Resonance  effects  will  not  be  considered,  nor  will  the  effects 
due  to  solar  pressure  or  third  bodies. 


Although  the  derivations  in  this  section  already  exist  in  the  literature,  in  the 


interest  of  completeness  they  are  presented  in  this  chapter. 


Derivation  q£  the  Geopotential  Equation 

Sir  Isaac  Newton  showed  that  in  inertial  space  the  gravitational  force  of  attraction 
between  two  bodies  can  be  written  as: 


F_ 

m 


(2.1) 


Newton  also  demonstrated  that  for  a  spherical  body  with  a  homogenous  distribution 
of  mass,  the  entire  mass  of  the  primary  body  acts  as  if  its  mass  existed  as  a  point 
particle  located  at  the  center  of  its  sphere.  If  a  planet  is  not  perfectly  spherical,  and/or 
does  not  have  a  homogenous  distribution  of  mass,  then  these  irregularities  will  effect  the 
motion  of  satellite  about  that  planet.  The  acceleration  which  a  satellite  experiences  (due 
to  the  mass  of  the  primary  body)  can  be  written  as  (19:49): 

a,  =  - 7 1' ( ,v ,  y,  z\  (2.2) 

The  i'i .v.y . ^  term  in  equation  (2.2)  can  be  solved  using  a  special  form  of  Poisson's 
Equation  that  is  known  as  Laplace’s  Equation  (the  derivation  of  these  equations  is  found 
in  Appendix  A)  which  in  cartesian  coordinates  is: 

72l  =  0  (2.3) 

The  equations  of  motion  of  a  satellite  in  orbit  around  a  planet  are  simpler  if 
expressed  in  spherical  polar  coordinates,  hence,  equation  (2.3)  becomes  (10:3): 


(2.4) 


/.  s. 


\_d_(  i  *(  ary _ i  a2r 

r2dr\  dr  )  r2  cos<pd<P\  ^  dp  )  r2cos2<pd\2 

where  r  =  radial  distance  from  the  center  of  the  attracting  body  to  the  satellite 


u  =  the  latitude 

t  =  the  longitude 

Equation  (2.4)  is  a  linear  partial  differential  equation  whose  solution  takes  the  form 

of  (10:4): 

1  r  ,  0  .  A.  =  R  r  <P  <p  A  K  (2.5) 

Because  i  r.o.\  describes  a  smooth  sphere  certain  boundary  conditions  must  be 
imposed  upon  equation  (2.5).  First,  in  order  to  prevent  a  jump  discontinuity  in  the  i 
function  it  must  have  the  same  value  at  .roi  and  ■>  2 n  .  Second,  to  prevent 
discontinuities  at  the  poles  of  the  sphere,  the  first  derivative  of  o  with  respect  to  ®  must 
equal  zero  when  ever  the  latitude  is  equal  to  odd  multiples  of  n/2. 


Substitution  of  equation  (2.5)  into  equation  (2.4)  yields: 


1  d  (  ,  d 

-  I  r  ‘  -  Rr<Ptp  1  A.  i  !  ,  , _ ^ 
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(2.6) 
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r 2  cos2  0  d \ ' 


R:n0[0,.lA  =  0 


The  right  hand  side  (RHS)  of  the  equation  (2.6)  can  be  written  in  terms  of  »  alone  as: 


c  o  s 2  0  elf  2dR 
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(2.7) 


Since  the  LHS  and  the  RHS  of  equation  (2.7)  are  independent,  set  them  equal  to 
the  constant  ‘ .  Hence  equation  (2.7)  implies: 


The  solution  to  equation  (2.8)  has  the  form: 


.  1  =  C  cos^ ;  A , 1 /2  A  +  Ssin  i  A' i 1/2  A.]  (--9) 

where,  in  addition  to  k  being  a  constant,  c  and  s  are  also  constants.  Further,  unlike  the 
constants  c  and  s,  k  can  not  be  an  arbitrary  value.  The  first  boundary  condition  in 
equation  (2.5)  implies  that  .*  must  be  equal  to  a  positive  integer.  Let  this  integer  be 
Hence  the  general  solution  to  equation  (2.8)  has  the  form: 

.lm  A  =  CmcosmA  +  Smsin  mA  (2.10) 

To  obtain  the  next  expression,  set  the  LHS  of  equation  (2.7)  equal  to  *  yielding: 


\  cl  (  2  cl R']  m2  1  df  cl<P\  (2.11) 

Rdr\  dr)  cos 2  <t>  cos  <t>  d  <p\  dip) 

Again  the  LHS  and  the  RHS  are  independent  of  each  other;  therefore,  equate  both  sides 
to  a  constant,  denoted  by  r.  Hence: 


l±(r*<** 

Rdr\  dr 


=  T 


(2.12) 


which  implies 
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(2.13) 


Here  the  second  boundary  condition  in  equation  (2.5)  imposes  the  already  mentioned 
conditional  values  of  /<•/#. 

Now  let  which  implies  - roso.r*.  From  this  the  mathematical  operator: 


is  derived.  Applying  this  operator  to  equation  (2.13)  yields: 
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2  d2<t>  J  m2  \  n 

cos 2  6 - -  -  <#>  — — —  7=0 
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(2.15) 


Therefore,  equation  (2.15)  becomes: 


-  -v 2  • - -  -  <t>  . - -  -7=0 

dx2  l  i  1  -  .v2 ,  J 


(2.16) 


Equation  (2.16)  is  the  algebraic  form  known  in  the  literature  as  Ferrier’s  form  of 


Legendre’s  Associated  Equation  whose  solution  has  the  form  (1:160-162): 


p?ix-.  =  i\  -x2r/2^pLix) 
dx 


(2.17) 


When  x-smo  equation  (2.17)  becomes: 


m 

P?  i  sin  0  '  =  !  cos  2  <t> ; 171/2 - P  L\ sin  <t> ) 

dxm 


(2.18) 


where  (1:132) 


„  ,  1  dL  r  2  L i 

P  L  \  s  i  n  <f> )  =  — - -  hsin  0-1 

2L  L  !  dlsin  0  J 


(2.19) 


In  equations  (2.17)  and  (2.18),  m  is  any  non  negative  integer,  and  t  is  an  integer  value 
which  is  the  number  of  times  p,isin«:  passes  through  zero  as  *  varies  from  0  to  n  (4:57). 

Equations  (2.18)  and  (2.19)  are  Rodrigues  formulas  giving  a  representation  of  the 
Legendre  polynomials.  An  alternate  expression  for  the  Legendre  polynomials  is  (1:132): 


therefore: 


Combining  equations  (2.10),  (2.21),  and  (2.28)  into  equation  (2.5)  yields  the  desired 
solution  to  equation  (2  4).  The  results,  equation  (2.29),  is  the  objective  of  this  section 
(19:55). 


r  ,  <t> .  A  =  - 
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(2.29) 


where  R,  =  the  equatorial  radius  of  the  primary  body 
g  =  the  universal  gravity  constant 

r  =  the  distance  from  the  center  of  the  primary  body  to  the  satellite 

The  c.„  and  the  terms  in  the  above  equation  are  constants  that  describe  the 
distribution  of  the  primary  body’s  mass,  and  are  known  in  the  literature  as  the  primary 
body’s  "gravity  model".  These  terms  are  dimensionless  since  the  dimensional  units  are 
carried  by  the  term  c.w. 

As  an  example  of  the  effects  of  a  primary  body’s  shape  and  distribution  of  mass 
upon  its  gravity  field,  an  18  by  18  gravity  model  for  the  planet  Mars  was  input  into 

equation  (2.29).  The  result  was  solved  for  values  of  latitude  and  longitude  that 

encompass  the  planet  at  an  altitude  of  500  km.  (see  the  program  Marsl  in  Appendix  E) 
The  results  are  plotted  in  Figure  2.1.  Note  the  checkerboard  or  "tesseral"  pattern  of 
alternating  regions  of  higher  and  lower  geopotential  than  would  exist  if  Mars  were  a 
perfect,  homogenous  sphere. 


The  Geopotential  Equation  as  a  Function  q£  ih£  Classical  Orbital  Elements 

Equation  (2.29)  is  given  in  spherical  polar  coordinates.  Since  satellite  motion  is 
often  described  in  terms  of  the  classical  orbital  elements,  it  is  necessary  to  derive 
equation  (2.29)  into  a  function  of  the  classical  orbital  elements.  This  derivation  is 
structured  on  the  work  of  Kaula,  Born,  and  Hildebrand,  see  references  10  and  3. 

Using  equations  (2.18)  and  (2.20)  rewrite  the  er(sm#>  term  in  the  above  equation  to 

yie'd: 
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(2.30) 


where 


P L  sin  0  = 
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4  2Lt'[L- tdi  L-2t\' 


Combining  equations  (2.30)  and  (2.31): 
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(2.32) 


Noting  that  ( 1:61 ) 
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Then 
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(2.34) 


Substituting  equation  (2.34)  into  equation  (2.32),  results  in: 


P 


m 
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(2.35) 


The  upper  limit  of  the  summation  changed  from  \l/ 2  due  to  the  denominator  term, 
t  ■  .  This  term  causes  any  value  of  1  greater  than  L-m  /2  to  make  the  factorial 
negative,  driving  the  factorial  to  infinity,  and  thus  driving  the  summation  to 

zero. 

Let  (10:6): 


1_^  2  L-  2  M 
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(2.36) 


Then  inserting  equations  (2.35)  and  (2.36)  into  equation  (2.29)  yields: 
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(2.37) 


In  order  to  utilize  Lagrange’s  Planetary  Equations,  equation  (2.37)  must  be 
rewritten  in  terms  of  the  six  orbital  elements,  a,  e,  1,  w,  n,  and  w.  Figure  2.2  shows  the 
relationships  between  the  various  angles. 


The  >  term  will  now  be  converted  into  the  orbital  elements.  From  Figure  2.2  it  is 
easily  seen  that  x-a-e;  however,  neither  a  nor  0  are  members  of  the  orbital  element  set. 
Therefore,  express  »  as: 


A  =  a-Q  -  0-0  =  a- O'*  0-0:  (2.38) 

where  *  =  the  longitude  of  the  projection  of  the  secondary  body  onto  the 

primary  body 


the  angle  from  the  x  axis  of  the  inertial  frame  to  the  longitude  of 
the  projection  of  the  secondary  body  onto  the  primary  body 

t>  =  the  angle  from  the  x  axis  of  the  inertial  frame  to  the  prime 


meridian  of  the  primary  body  (also  known  as  the  "local  sidereal 
time") 


Figure  2.2 


Applying  equation  (2.38)  to  the  cos mk  and  s.nm*  terms  of  equation  (2.37)  yields: 

cor,  m  K  =  cos  m  a  ~  Q  *  m  n  -  0  (2.39) 

sm  rn  A  -sin  m  a-Q  *  Q-9  (2.39a) 


Applying  equations  (Bl.l)  and  (B1.2),  found  in  Appendix  B  to  equations  (2.39)  and 
1 2.39a)  yields: 
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( 2  40a ) 


Noting  the  angular  relationships  in  Figure  2.2,  and  using  the  properties  of  spherical 
trigonometry  the  following  relationships  are  evident: 


cos  a  -  O  = 
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(2.41) 


sin  a  -  O  =  tan  0cot i 


sm  0  =  sin  co  *  /  'sin  i 


(2.42) 


(2.43) 


Looking  at  the  .-os  a- a  terms  of  equation  (2.40)  and  applying  equations  (B2.5),  (B2.7). 


(2  41 ),  and  (2.42)  yields: 
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I  he  same  process  applied  to  the  sine  terms  of  equations  (2.40)  yields: 
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Injecting  equations  (2.44)  and  (2.45)  into  equations  (2.40)  will  result  in: 


(2.44) 


„  _  v-  (  rn\  cos’"  ’  co*  /  sin’  to*  /  cos  ’  i 

os  rn  A  R  h  >  \)  -  -  -  - 

fr0{s  r  cos'"1  0 

x  cosm  0-0  *  ysm  m  O  -  0  ] 


i  rvamft  M 


i 


m3  cos"'  '  (i)  *  /  ‘.in'  a>  *■  f  cos'i 
;i  n  rn  A  =  [  RI ■  ■  .  i  y '  ,„ 

“V‘-'  cos  0 


(2.47) 


x  smm  0-0  -  yrosm  0-0  ] 

Applying  equation  (2.36),  equation  (2.35)  can  be  written  as: 
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Inject  equation  (2.43)  into  equation  (2.48)  to  yield: 
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Substituting  equations  (2.46),  (2.47),  and  (2.49)  into  equation  (2.29)  yields: 
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The  last  term  of  equation  (2.50)  is  in  the  form  of  smscoss.  Applying  equations  (B3.2) 
and  (B3.3),  from  Appendix  B,  to  equation  (2.50)  yields: 
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terms  of  equation  (2.51)  yields  (3:5): 
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Since  ■  in  equation  (2.52)  is  a  real  physical  quantity  it  is  necessary  to  determine  the  real 
part  of  the  bracketed  term.  Consider  the  r  term. 


jS  ^  j  l  rn  2 1  *  s  _^s  |  /  j  i  m  2  ( •  i  _  ^  s  ^  L  •  m  •  ?  t  s 


(2.53) 
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Remember  that  *  is  the  integer  part  of  L-m  /2%  so  if  r-m  is  odd  (3:7): 
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2  2 


(2.54) 
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When  /  -  ™  is  even: 


L  -  rn 


(2.56) 


j%  -  j  L  m " 21  *  * 


i  ‘m=  -i 


(2.57) 


As  a  result  of  equations  (2.55)  and  (2.57)  equation  (2.52)  becomes: 
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(2.58) 


->  im  l  -  m . even 
C  1  -  in  .  od  d 


sin  L-2t-2c-2d  to  *  f  +  m  0-0  > 


Transform  equation  (2.58)  so  terms  of  the  form: 


l  -  2  p  uj  +  f+m  0-0 


I 


can  be  collected  together.  This  is  accomplished  by  letting  (3:8): 


p  m  t  *c  *  d 


(2.59) 


This  implies 


/.  -  2t  -  2.c  -  2d  =  L  -  2  p 


(2.60) 
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Hence,  equation  (2.58)  becomes: 


.  •  0  m  •  0  \  p  /  f  -  0 
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(2.61) 
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Evaluating  equation  (2.36)  yields  the  following  relationship: 


0<t<k  (2.62) 

Likewise,  an  evaluation  of  the  binomial  coefficient  terms  of  equations  (2.58*)  and  (2.61) 

>  ields: 

0  £  s  <  m  (2.62a) 

O  <  c  <  L  -  m  -  2t  *  s  (2.62b) 

0  <  d  <  m  -  s  (2.62c) 

0<p<L  (2.62d) 

However,  according  to  equation  (2.59)  i-p-c-d.  Since  both  c  and  d  have  minimum 
values  of  zero,  t„..  is  equal  to  p.  This  implies  that  the  maximum  value  of  t  will  be  the 
smaller  value  of  p  or  t,  and  equation  (2.62)  becomes: 

o <t<  the  smaller  of  *  or  p  (2.62e) 


Grouping  selected  terms  from  equation  (2.61),  and  taking  into  account  the  possible 
values  for  \  ,  r.  and  p  from  equations  (2.62),  leads  to  the  definition  (10:34): 
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l.  -  m  -  It  *■  s\f  m  -  s 
p  -  (  -  c 


is  known  in  the  literature  as  the  Inclination  Function.  A  table  of  values  for  this 
function  is  given  in  Appendix  C.  Rewriting  equation  (2.61)  in  terms  of  the  inclination 
function  gives: 
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(2.64) 


5 (<n  L  -  m  .pi  en 
C Lm  L  -  m  .  odd 


sin[(Z-2p)(tu  +  /  )+m(il-0)] 


Where  equation  (2.64)  is  in  a  form  that  is  for  a  particular  value  of  L  and  m. 

Next,  equation  (2.64)  must  be  written  so  that  the  r  and  /  terms  are  expressed  in 
terms  of  i,  v,  and  a.  From  equation  (2.64)  isolate  any  particular 


1  cos 


Z  -  2  p  cu*f  *mQ-Qi 


(2.65) 


term  and  let 


f  =  L  -  2  p  )o>  +  m  fl  -  6 , 


(2.66) 


Equation  (2.65)  becomes: 


TTT  C°S  l.  -  2  p  f  *  e 
rL  ’  L  sin 


Now,  consider  the  term: 


(2.67) 


t  cos  Z-2p  f*  f  =  cos.  Z  -  2  p  : /  ’ cos e  -  sin[ ‘  Z  -  2  p  1  /  jsin  f 

Here  Born  et  al.  introduces  the  term: 


(2.68) 


^  j  exp  jmf  =  Y "  exp  iyimi 
Where  x  '  is  known  as  Hansen’s  coefficients  (2:2): 
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in  the  above  equations  is  the  Eccentric  Anomaly.  Employing  equation  (2.69)  along 
with  the  following  relationships: 


o\p  jmf  =  co smf  +  ]  sin  mf 


o\p  jim  =  cos  im  *  jsmim 


i  =!  -  7  p  *  q  m  =  l  -  2  p  rt  =  -  L 


and  noting  that  from  equation  (2.69)  follows: 
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(2.73) 
(2.73a) 

(2.74) 

(2.75) 

(2.75a) 


changes  equation  (2.65)  into: 
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(2.76) 
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L  -  2  p  oj  *  L  -  2  p  +  q  \f  +  m  O  -  Q 


Next  a  determination  of  the  characteristics  of  Hansen's  coefficients  is  necessary.  It 
is  beyond  the  scope  of  this  chapter,  but  it  can  be  shown  that  for  the  case  at  hand  (2:5): 
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then  direct  substitution  into  equation  (2.77)  yields: 
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If  ■'  -  .  ,  then  equation  (2.77)  yields: 
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Here  Born  et  al.  makes  the  following  definition: 


p  =  p  for  p  <  L/2 


p'  =  L  -  p  for  p>  L/2 


This  implies: 
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Substituting  directly  into  equations  (2.83)  and  (2.84)  yields: 
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By  examining  the  combination  of  cases  for  ?>o,  <j<o,  p<l/ 2,  and  p>l/2  it  can  be  shown 
that  (3:14-15): 


X  -  L  1  .  L 

A  L  2p-  < 


2p 


'  Z.p<7 


P  =  !  -  1 


P 2  Lp 


zL  ^  LPQt  Q  Lpqk  P 


(2.88) 


where 
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In  equations  (2.89)  and  (2.90)  the  following  conditions  hold,  n-k-q  if  q  >  o.  If  q  <o 
then  h -k.  Also,  p  ~p  and  q  -q  if  p<l/2.  If  p>L/ 2  then  p  -l~p  and  q--q. 

The  e  1  term  is  known  as  the  Eccentricity  Function.  A  list  of  this  function's 
values  is  in  Appendix  C. 

Equations  (2.63),  (2.64),  and  (2.88)  allow  equation  (2.29)  to  be  written  into  the 

form: 
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where 
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A  satellite  moving  through  an  atmosphere  experiences  a  force  perpendicular  to  its 
flight  path  ("lift"),  and  a  force  in  the  opposite  direction  to  its  flight  path  ("drag"). 
Because  of  variations  in  a  satellite’s  attitude,  the  resultant  lift  force  is  usually  zero.  This 
is  especially  true  for  spherical  satellites,  or  satellites  whose  length  is  greater  than  its 
diameter.  Even  if  the  resultant  lift  force  is  not  zero,  its  effects,  when  compared  with 
drag,  are  still  negligible  (15:295).  Drag,  on  the  other  hand  can  have  a  profound  effect 
on  the  orbit  of  3  satellite. 

In  this  analysis,  the  atmosphere  is  modeled  as  a  locally  exponential  atmosphere. 
Therefore,  the  density  of  the  atmosphere  is  decreasing  exponentially  with  altitude, 
implying  that  drag’s  predominant  effects  occur  when  the  satellite  is  near  its  closet 
approach  to  a  planet.  At  this  point  the  flight  path  angle  is  approximately  zero.  Thus, 
drag  will  be  acting  directly  opposite  to  the  satellite’s  velocity  vector.  This  will  have  the 
effect  of  slowing  the  satellite  down,  and  hence,  decreasing  its  energy.  The  decrease  in 
the  satellite’s  energy  will  result  in  a  decrease  in  the  semi  major  axis  a,  and  the 
eccentricity,  *.  Although  periapsis  altitude  will  decrease  somewhat,  this  decrease  is  very 
small  when  compared  with  the  resulting  decrease  in  the  apoapsis  altitude.  The  over  all 
effect  of  drag  will  to  be  to  "circularize"  the  orbit. 

If  the  atmosphere  were  perfectly  spherical  and  nonrotating,  the  reduction  in  a  and 
would  be  drag’s  only  effects  on  the  orbit.  However,  atmospheres  share  the  same 
propensity  for  oblateness  as  their  planets  and  tend  to  rotate.  The  oblateness  of  the 
atmosphere  will  induce  small  changes  in  the  argument  of  periapsis,  u,,  while  the  rotation 
of  the  atmosphere  results  in  small  lateral  forces  on  the  satellite.  These  lateral  forces 
cause  increasing  changes  in  the  angle  of  inclination,  i,  and  small  periodic  changes  in  the 
longitude  of  the  ascending  node,  n  (11:6-7). 


-s  *> 


The  atmospheric  drag  on  a  satellite  may  be  expressed  as  (15:295): 


\ 


D  = 


pl2SC0 

2m 


where  D 

p 

v 
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the  force  of  drag 
the  atmospheric  density 

the  velocity  of  the  satellite  relative  to  the  atmosphere 
the  effective  area  of  the  satellite 
the  coefficient  of  drag 
the  mass  of  the  satellite 


(3.1) 


Atmospheric  Density 

Appendix  D  develops  the  expression  used  for  a  locally  exponential  atmosphere. 
This  expression  is: 


P - P oe*P 


9mA 

RT^j 


(3.2) 


fn  equation  (3.2),  is  equal  to  the  altitude  above  the  planet.  To  write  equation 
(3.2)  in  terms  of  the  radial  distance  (r)  from  the  center  of  the  planet  let: 

T-  -  r  -  R  p  (3.3) 

where  «,  is  the  radius  of  the  planet.  Using  an  expression  for  the  rectangular 
components  of  a  point  on  the  surface  of  a  planet  as  found  in  Escobal,  page  26,  r,  can  be 


written  as: 


where 


the  equatorial  radius  of  the  planet 


.  =  the  eccentricity  of  the  planet’s  shape 

»  =  the  latitude 

Applying  equation  (3.3)  to  equation  (3.2)  yields: 

P-P„e*r(-{af):r-S,  ) 


(3.3) 


The  bracketed  term  is  equa.  to  i/h,  where  h  is  the  "scale  height"  and  is  equal  to  the 
change  in  altitude  required  in  order  for  the  density  to  change  by  one  exponential.  As 
can  be  seen  from  equation  (3.5)  it  is  not  constant;  however,  at  the  altitudes  that  the 
satellite  will  experience  significant  air  drag  «  is  so  large  it  can  be  treated  as  a  constant. 
For  example,  using  data  obtain  from  the  Viking  1  space  craft,  at  200  km  altitude  the 
scale  height  is  14.1387  km  (16:4368-4373).  It  is  because  the  scale  height  can  be 
considered  a  constant  over  some  small  altitude  band  that  the  assumption  of  a  locally 
exponential  decreasing  atmosphere  may  be  made  (18:4).  Therefore,  equation  (3.5) 
becomes: 

/  r-R„  \  (3.6) 


Velocity  With  Respect  t£  ih£  Atmosphere 
Let. 

f  =  velocity  of  satellite  relative  to  the  atmosphere 
r  =  velocity  of  the  satellite  relative  to  the  planet 

u  =  velocity  of  the  atmosphere  relative  to  the  planet  (atmosphere  assumed  to 

be  moving  west  to  east) 


Figure  3.1  shows  the  angular  relationships  between  these  vectors. 


Angular  Relationships  Between  the  Different  Types  of  Velocities 

Figure  3.1. 


From  Figure  3.1  it  can  be  seen  that: 


l  “  v  -  u 


Applying  the  law  of  cosines  yields: 


1  2  =  i  2  *  u 2  -  2iu  cos  y 


Assume  that  the  atmosphere  rotates  with  an  angular  velocity  w>  about  the  planet.  Then 

u  =  r  oj  cos<p  (3.9) 

where  r  =  radial  distance  from  the  center  of  the  planet 
«  =  latitude 

Using  spherical  trigonometry  and  the  angular  relationships  in  Figure  3.1: 

cos  i  =  cos  0  cos  y  '  (3.10) 

This  analysis  assumes,  since  the  most  profound  effects  occur  at  periapsis,  that  the 
satellite  is  at  its  periapsis  point.  Thus  implying  that,  y  -y.  y  is  still  a  good 
approximation  for  y  even  when  the  satellite  is  not  at  its  periapsis  point.  However,  the 
satellite  must  be  with  in  two  scale  heights  of  periapsis  altitude  to  keep  the  error  of 
assuming  v  -  y  to  less  than  one  percent  (11:23). 

Therefore,  assuming  v  -y  and  applying  equation  (3.10)  to  equation  (3.9)  yields: 

u  =  r  to  cos  0  (3.11) 

cos  i 

=  rcu 

cos  y 

lie  os  y  =  r  oj  c  os  t 

Substituting  equations  (3.9)  and  (3.1  1 )  into  equation  (3.8)  yields: 

l  r~uj  V  ,  2  ,  (3.12) 

1  =  (  I  I  -  —  cos  t  j  -  r  ‘  cu  cos  0  -  cos’; 

For  the  planet  Mars,  the  atmosphere  rotates  with  approximately  the  same  angular 
velocity  as  the  planet  ( 18.3).  Therefore,  *>-  oooososo  radians  per  second  (13:2-3).  This 
small  value  for  <■  results  in  the  r< term  in  the  3bove  equation  being  vanishingly  small 


when  compared  to  1  and  will  be  neglected.  Further,  since  drag  effects  the  periapsis 


altitude,  velocity,  and  angle  of  inclination,  equation  (3.12)  must  be  rewritten  for  some 
reference  periapsis  altitude,  velocity,  and  inclination.  This  is  accomplished  by  letting 
,  >-t.  .and  1-1  Equation  (3. 12)  becomes: 
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(3.13) 


The  Cross  Sectional  Area.  S 


The  cross  sectional  area  effecting  drag,  S,  will  be  a  function  of  the  satellite’s  shape 
and  flight  path  angle.  Due  to  the  array  of  scientific  sensors  desired  for  a  Mars  mission, 
the  satellite’s  shape  will  most  likely  be  very  irregular,  implying  that  the  effective  cross 
sectional  area  may  not  be  known.  No  matter  what  the  shape,  a  satellite  in  uncontrolled 
flight  will  have  a  tendency  to  rotate  about  its  axis  of  maximum  moment  of  inertia 
( 8:369- 371).  For  cylindrical  shaped  satellites  with  a  length  to  diameter  (L/d)  ratio 
greater  than  roughly  2,  this  rotation  will  cause  the  satellite  to  move  through  the 
atmosphere  tumbling  end  over  end,  or  revolving  like  an  aircraft  propeller  (11:16). 

In  the  first  case  a  mean  value  of  s  is: 


(3.14) 


and  in  the  second  case: 


5  =  Id  (3.15) 

where  /  =  the  length  of  the  satellite 

t  =  the  diameter  of  the  satellite 

If  the  direction  of  the  spin  axis  is  not  known,  then  the  mean  value  of  s  is 
somewhere  in-between  the  values  given  in  equation  (3.14)  and  (3.15).  Averaging  these 
two  equations  yields: 
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This  \alue  will  never  be  more  than  15  percent  off  the  extreme  case  (satellite  spinning 
like  a  propeller). 


When  L  d  is  less  than  1/2,  the  spin  axis  becomes  the  axis  of  symmetry  In  this 
case,  if  the  spin  axis  is  aligned  with  the  satellite’s  direction  of  motion: 

S  =  rr  r 2  <3.17) 

If  the  spin  axis  is  perpendicular  to  the  flight  path,  :>  is  given  by  equation  (3.15).  In  this 
case,  if  ;  is  much  smaller  than  •<  the  value  of  s  can  become  very  small.  This  implies  that 
the  error  associated  in  averaging  the  values  of  equations  (3.15)  and  (3.17),  when  the 
direction  of  the  spin  axis  is  unknown,  can  yield  differences  between  the  actual  and 
estimated  values  of  s  that  are  much  greater  than  those  of  the  previous  case. 

For  this  thesis,  based  on  a  rough  estimate  on  the  size  of  satellites  currently  orbiting 
the  earth,  a  cross  sectional  area  of  s-  io-n-  will  be  used. 

The  Coefficient  qf  Drag.  Cp 

The  coefficient  of  drag  is  dependent  upon  the  density  of  the  atmosphere,  the 
Reynolds  number,  angle  of  attack,  the  shape,  and  the  speed  of  the  satellite.  These 
parameters  not  only  vary  from  satellite  configuration  to  satellite  configuration,  but  can 
also  vary  through  out  the  satellite’s  flight  path. 

As  the  density  increases  three  distinct  regions  of  atmospheric  flow  are  encountered. 
First,  continuum  flow,  is  the  region  where  the  atmosphere  deforms  continuously  under 
the  shear  force  applied  by  the  moving  satellite.  The  Viking  project  found  that  for  Mars 
this  region  exist  from  the  surface  to  about  90  km  altitude.  For  Viking  the  coefficient  of 
drag  in  this  region  was  approximately  1.47.  Next,  the  slip  flow  region,  which  exist  from 
about  90  km  to  115  km,  is  a  region  of  transition  between  continuum  flow  and  free 
molecular  flow,  the  third  region.  Free  molecular  flow  exist  when  the  distance  that  a 


molecule  can  travel  with  out  striking  another  molecule,  its  mean  free  path,  i  greater 
than  the  dimensions  of  the  satellite.  For  Mars  this  region  exist  for  altitudes  greater  than 

r-'UghK  115km 

This  thesis  is  concerned  primarily  with  the  region  of  free  molecular  flow.  In  this 
region  the  coefficient  of  drag  can  \ary  as  the  angle  of  attack  of  the  satellite  varies,  and 
will  he  on  the  order  of  2  0  to  2.25.  A  Co  of  2.0  will  be  used  in  this  thesis.  This  value 
was  chosen  because  it  was  the  coefficient  of  drag  used  on  the  Viking  mission  (16:4369). 


IV.  Computer  Program  Validation 


Description  of  the  Program 

Part  of  the  analysis  of  this  thesis  was  carried  out  using  the  Artificial  Satellite 
Analysis  Program  (ASAP),  see  reference  13.  This  program  uses  Cowell’s  method. 
Essentially  this  involves  taking  the  state  vector  of  the  satellite  with  respect  to  an  x,  y,  z 
coordinate  system  whose  origin  is  at  the  center  of  the  central  body,  whose  xy  plane  lies 
in  the  plane  of  the  equator,  and  whose  z  axis  goes  through  the  north  pole  of  the  central 
body;  and  then  solving  the  associated  equations  of  motion  via  a  numerical  integration 
package.  ASAP  uses  an  8th  order  Runge-Kutta  integrator  that  requires  the  equations  of 
motion  be  written  as  a  set  of  first  order  differential  equations.  This  process  looks  like 
i  13:3-1); 

$  =  v  .  y  .  \  ,  y  .  r.  r  (4.1) 

where  ■  =  i  ,  =  velocity  in  the  X  direction 

-  =  i  .  =  velocity  in  Y  direction 

=  >  ,  =  velocity  in  Z  direction 

Applying  Newton's  second  law  (to  determine  the  equations  of  motion),  and  keeping  in 
mind  that  the  Runge-Kutta  package  used  requires  a  set  of  first  order  differential 
equations  yields  (13:3-1): 

,  •  -v  „  (4.2) 
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where  u  =  the  universal  gravity  constant  multiplied  by  the  mass  of  the  central 
body 


This  thesis  will  only  consider  those  perturbing  effects  caused  by  the  central  body  and 
atmospheric  drag. 

Perturbations  due  &  th£  Central  Body.  ASAP  uses  equations  (2.29),  (2.91),  and 
(2.92)  in  order  to  find  the  geopotential  in  terms  of  latitude,  longitude,  radial  position, 
and  the  classical  orbital  elements,  i,  e,  a.,  w,  and  n.  The  implementation  of  these 
equations  into  a  form  acceptable  to  the  Runge-K.utta  integrator  requires  the  conversion 
to  cartesian  coordinates.  This  is  accomplished  by  rewriting  equation  (2.29)  such  that 
only  effects  due  to  departures  in  the  central  body’s  shape  from  a  perfect,  homogenous 
sphere  are  considered.  The  resulting  equation  is: 
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The  perturbing  portion  of  equations  (4.2)  through  (4.4)  due  to  the  central  body  can  now 
be  w  ritten  as  ( 1 3:3- 1 ): 
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Perturbations  due  to  Atmospheric  Drag.  Since  atmospheric  drag  acts  to  retard  the 
motion  of  a  satellite,  the  equation  of  motion  of  atmospheric  drag  used  by  ASAP  is  the 
negative  of  equation  (3.1).  As  a  model  for  the  atmospheric  density  change  with  altitude, 
equation  (3.6)  is  employed;  however,  the  selection  of  a  reference  height  from  which  to 
base  density  calculations  is  allowed.  This  is  implemented  by  replacing  the  Rr  term  in 
equation  (3.6)  with  a  reference  height  term,  h:.  The  program  also  takes  into  account  the 
departure  in  the  central  body's  shape  from  a  perfect  sphere  by  use  of  equation  (3.4). 


Program  Validation. 

The  following  equations  were  used  in  the  validation  of  ASAP.  They  were  derived 
from  the  Lagrange  Planetary  equations  where  the  disturbing  function  is  derived  from 
equation  (2.91),  using  only  the  zonal  harmonics  up  to  and  including  the  6th  zonal 
(14:28-30): 
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Each  of  the  above  4  orbital  elements  (<>,  i,  n,  u>)  were  analyzed  and  predictions  were 
made  as  to  what  values  will  drive  the  change  in  each  element  to  zero.  These  predictions 
were  then  tested  by  running  ASAP  with  the  appropriate  elements.  If  ASAP  is  reliable 


( 


both  the  predictions  and  the  ASAP  output  should  agree. 

Because  ^  contributes  only  to  the  long  term  perturbations,  which  are  periodic  over 
one  axial  period  (the  time  for  the  line  of  apsides  to  make  one  complete  revolution),  all 
trigonometric  terms  containing  cu  will  be  set  equal  to  zero.  This  greatly  simplifies  the 
above  equations,  and  is  valid  due  to  the  method  of  averaging  when  applied  to  the  ^ 
terms  of  equation  (2.91).  Setting  <*>  equal  to  zero  causes  all  the  odd  zonal  harmonics  in 
equations  (4.12)  through  (4.39)  to  go  to  zero,  and  eliminates  many  other  terms  from  the 
even  zonal  harmonics. 

Eccentricity  and  Inclination.  Equations  (4.12)  through  (4.25)  do  not  have  any  non 
zero  terms  once  trigonometric  functions  of  ^  have  been  set  to  zero.  This  implies  that  no 
matter  what  the  size,  shape,  or  orientation  of  the  orbit  the  secular  changes  in 
eccentricity  and  inclination  due  to  zonal  harmonics  are  zero.  Eccentricity  and 
inclination  will  experience  a  short  term  change  due  to  a  change  in  the  mean  anomaly, 
and  also  a  long  term  change  due  to  precession  of  the  line  of  apsides;  however,  since 
both  these  effects  are  periodic,  and  since  the  change  in  w>  over  one  orbital  period  is 
small  compared  to  the  change  in  mean  anomaly,  the  change  in  eccentricity  and 
inclination  over  one  orbital  period  will  be  almost  zero  while  the  change  in  eccentricity 
and  inclination  over  on  axial  period  will  be  zero.  This  prediction  is  also  supported  by 
Roy.  page  290. 

Several  computer  runs  were  made  with  ASAP  using  different  input  values.  These 
data  runs  considered  the  perturbative  effects  due  to  zonal  harmonics  up  to  and  including 
an  order  of  six.  In  all  cases  the  output  was  consistent  with  the  above  predictions. 
Figures  4.1  through  4.4  are  a  representative  sample  of  the  output,  and  indicates  the 
change  in  eccentricity  and  inclination  over  one  orbital  period,  and  one  axial  period. 
Table  4.1  lists  the  input  orbital  elements  used  to  generate  Figures  4.1  through  4.10. 


*.  A 

«VV 


/•  /.  «•_  ** .  y* 


S2.260C 


! 


Change  in  Inclination  Over  One  Axial  Period  (6X0  Gravity  Field) 

Figure  4.4 

Note,  due  to  inaccuracies  in  calculating  the  exact  axial  period,  the  change  in  eccentricity 
and  inclination  shown  in  Figures  4.3  and  4.4  are  not  exactly  zero. 

Longitude  of  the  Ascending  Node.  Equations  (4.26)  and  (4.32)  have  the  term  cosi 
as  a  common  denominator.  Thus,  any  value  of  the  inclination  that  drives  the  cosi  term 
to  zero  will  cause  the  change  in  the  Longitude  of  the  Ascending  Node  (n)  to  also  equal 

zero.  A  polar  orbit  (<-9o)  has  long  been  known  to  yield  an  equal  zero.  Figure  4.5  shows 

the  ASAP  output  given  an  input  value  of  i  equal  to  ninety  degrees.  As  expected, 

throughout  the  orbital  period  there  is  no  change  in  n. 
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Change  in  Ascending  Node  for  a  Polar  Orbit  (6X0  Gravity  Field) 

Figure  4.5 

In  addition  to  the  predominant  cos<  term,  equations  (4.27)  through  (4.32)  also 
contain  other  trigonometric  functions  of  i.  A  search  for  values  of  <  (other  than  i  -<?o, 
(-  270)  that  will  cause  jo  to  equal  zero  was  made  by  inputting  equations  (4.27)  through 
(4.32)  into  equation  (4.26).  The  cosi  terms  were  eliminated  by  setting  jo  equal  to  zero 
and  then  dividing  by  cosi.  The  only  i  terms  left  in  the  equation  are  powers  of  sin. 
Terms  were  grouped  by  the  power  of  their  associated  sim  terms  thus  producing  a  4th 
degree  polynomial.  By  making  the  change  of  variable  the  polynomial  is  reduced  to 

a  quadratic.  This  quadratic  was  solved  using  the  computer  program  Capmega  given  in 
Appendix  E.  The  results  show  that  for  eccentricities  from  0  to  .9,  and  for  inclinations 
from  0  to  90,  there  is  no  other  value  of  <  that  will  yield  jo  equal  to  zero  other  than 
those  values  of  <  associated  with  a  polar  orbit. 

Note  that  equations  (4.21)  through  (4.25)  also  have  a  common  denominator  of  cosi, 
thus  implying  that  the  change  in  inclination  will  also  be  zero  if  in  a  polar  orbit.  This 
gives  another  opportunity  to  validate  ASAP  by  noting  its  predicted  change  in  inclination 
for  a  polar  orbit.  Figure  4.6  shows  the  results  of  this  procedure. 
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Change  in  Inclination  for  a  Polar  Orbit  (6X0  Gravity  Field) 

Figure  4.6 

Argument  of  the  Periaosis.  In  a  process  similar  to  the  one  described  above, 
equations  (4.34)  through  (4.39)  were  substituted  into  equation  (4.3 3),  and  terms  of 
similar  powers  of  smi  were  grouped  together.  The  program  Omega  (found  in  Appendix 
E)  was  used  to  solve  the  resulting  polynomial.  The  results  yielded  a  particular  value  for 
,  taking  into  account  the  zonal  harmonics  up  to  order  six,  that  causes  Jo>  to  equal  zero. 
This  value  is  know  as  the  critical  value  of  i,  and  is  dependent  on  the  semi  major  axis 
and  the  eccentricity  of  the  orbit.  A  critical  value  of  <  was  determined  for  several 
different  values  of  eccentricity  and  semi  major  axis.  These  values  were  inputted  into 
ASAP.  In  each  case  ASAP  yielded  the  proper  result. 

As  with  the  eccentricity  and  inclination,  <*>  is  subject  to  short  and  long  term 
perturbations;  therefore,  when  the  inclination  is  at  its  critical  value,  the  change  in 
o%er  one  orbital  period  should  be  close  to  zero,  while  the  change  over  one  axial  period 
should  be  exactly  zero.  Figure  4.7  shows  that  the  change  in  ^  over  one  orbital  period  is 
indeed  almost  zero.  Figures  4.8  through  4.9  demonstrate  the  closed  nature  of  the  change 
in  eccentricity,  and  inclination  vs.  the  change  in  the  argument  of  the  periapsis  over  one 
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Figure  4.10 

A  look  at  the  change  in  ~u  over  one  axial  period  was  not  made  because  of  the 
extreme!)  long  axial  period  associated  with  the  critical  values  of  i  (on  the  order  of  15 

>  ears  > 

Atmospheric  Drag.  Equation  (4.40)  describes  the  change  in  the  semi  major  axis 
due  solely  to  atmospheric  drag  over  one  orbital  period  (11:41): 
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liquation  (4  44)  describes  the  change  in  eccentricity  due  solely  to  atmospheric  drag  over 
one  orbital  period  (11:41): 
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Putting  equation  (4.44)  into  integral  form  yields: 
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(4.45) 


Fquations  (4.40)  and  (4.45)  where  solved  using  an  8th  order  Gaussian-Legendre 
quadrature  method  (see  program  Dsemi  in  Appendix  E),  and  the  resulting  output 
compared  to  ASAP.  The  results  showed  that  the  above  equations  and  ASAP  give 
reasonablv  close  answers. 


V.  Analysis 

The  Mars  Geoscience  Climatology  Phasing  Orbit 

The  Mars  Geoscience  Climatology  Orbiter  (MGCO)  phasing  orbit  is  a  frozen  orbit 
planned  for  the  next  L'.S.  space  mission  to  Mars.  Table  5.1  list  the  elements  of  this 
orbit  The  Longitude  of  the  Ascending  Node  (o)  of  the  actual  orbiter  will  be  set  by  the 
approach  asymptote,  which,  for  the  purposes  of  this  analysis  will  be  90  degrees. 

Orbital  Elements  for  the  MGCO  Phasing  Orbit  (17:3) 

Table  5.1 


Input  Or¬ 
bital  Ele¬ 
ments  for: 

a  km 

P 

MGCO 

Phasing 

Orbit 

3747.2 

0.0081 

degrees 

90.00 


n  de¬ 
grees 

w  de¬ 
grees 

v  de¬ 
grees 

90.00 

270.00 

90.00 

In  order  to  determine  the  predominant  characteristics  of  a  frozen  orbit  the  above 
elements  were  inputted  into  ASAP,  propagated  for  one  and  three  orbital  periods,  and  for 
one  axial  period  using  both  a  6  X  0  and  a  6  X  6  gravity  field  (with  and  without 
ltmospheric  drag)  The  axial  period  was  estimated  by  using: 

<lu>  3  J  2R\  f  b  ,  \  P  (5  I) 

-II  7  a  l-A  •V'2S,n  0\  »> 


.  it :  ■  n  i5  I  i  was  derived  from  the  Lagrange  Planetary  Equations  using  only  the  second 
il  harnv>mc  of  equation  (2.91).  Analysis  on  the  output  revealed  that  over  one  orbital 
!  rhe  atmospheric  drag  has  no  appreciable  effect.  However,  for  a  6  X  0  gravity 
!  ’h-  Increase  in  the  semi  major  axis  over  one  axial  period  is  0.5  meters  more  when 
c  ;  <-nr  rhan  when  absent  For  a  6  X  6  gravity  field,  atmospheric  drag  causes  the 
i  r  avis  to  decrease  bv  0  1  meters  more  than  the  presence  of  a  6  X  6  gravity 
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Figures  5.1  through  5.3  show  the  effects  due  to  atmospheric  drag  and  a  6  X  0 
gravity  field,  while  Figures  5.4  through  5.6  show  the  effects  due  to  atmospheric  drag 
and  a  6  X  6  gravity  field.  Figures  5.1  and  5.2  are  'losed  curves,  asserting  that  the 
\alues  of  the  argument  of  the  periapsis,  the  eccentricity,  and  the  semi  major  axis  are 
bounded.  In  examining  Figure  5.3,  it  should  be  remembered  that  the  inclination  does 
not  change  over  one  orbital  period  for  polar  orbits  (see  equations  (4.19)  through  (4.25)); 
therefore.  Figure  5.3  shows  that  the  values  of  the  argument  of  the  periapsis  and  the 
angle  of  inclination  are  also  bounded.  This  bounded  condition  implies  that  the  values  of 
the  argument  of  the  periapsis,  the  eccentricity,  the  semi  major  axis,  and  the  inclination 
are  periodic  over  one  orbit.  This  situation  changes  when  a  6  X  6  gravity  field  is 
introduced.  Figure  5.4  reveals  that  the  initial  value  and  the  final  value  of  both  the 
argument  of  the  periapsis  and  the  eccentricity  are  not  the  same.  Over  one  orbital  period 
the  argument  of  the  periapsis  changes  from  to  =  275.67128  degrees  to  <o  =  276.55445 
degrees,  a  change  of  approximately  .32  percent  over  the  initial  value.  The  eccentricity 
changes  from  »  =  .00810551  to  e  =  .00806312,  representing  a  change  of  .5  percent  over 
the  initial  value.  Figure  5.5  reveals  that  the  the  semi  major  axis  changes  from  a  = 
3756.23351  km  to  a  =  3756.16521  km,  giving  a  change  of  approximately  .00182  percent. 
Although  Figure  5.6  shows  no  discernible  difference  from  Figure  5.3,  an  analysis  of  the 
data  shows  that  there  is  a  0.04756  degree  change  in  inclination  over  one  orbital  period 
when  in  a  6  X  6  gravity  field.  These  changes  in  the  orbital  parameters  indicate  that  the 
above  orbital  parameters  are  not  periodic  over  one  orbital  period  when  in  the  presence 
of  a  6  X  6  gravity  field.  To  test  these  conclusions,  the  MGCO  phasing  orbit  was 
propagated  over  three  orbital  periods  for  a  6  X  0  gravity  field  and  a  6  X  6  gravity  field, 
both  with  drag.  For  a  6  X  0  gravity  field,  Figures  5.7  through  5  9  show  that  the  orbit 
continues  to  exhibit  the  same  periodic  behavior  in  the  argument  of  the  periapsis,  the 
eccentricity,  the  semi  major  axis,  and  the  inclination  over  three  orbital  periods  as  was 
established  in  the  first  orbit.  This  confirms  the  predictions  made  from  Figures  5.1 
through  5.3. 
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Figure  5.4 
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In  the  next  step  ot'  this  analysis,  the  phasing  orbit  was  propagated  over  one  axial 
period.  Figure  5  13  shows  the  effect  of  a  6  X  0  gravity  field,  with  drag.  The  figure 
indicates  that  the  \alues  of  the  argument  of  the  periapsis,  and  the  eccentricity  are  not 
quite  periodic  over  the  axial  period,  which  is  not  correct.  The  axial  period  was 
estimated  from  Equation  (5.1),  which  does  not  take  into  account  zonal  harmonics  greater 
than  two,  nor  any  of  the  sectoral  harmonics;  therefore,  it  does  not  return  the  exact  axial 
period  If  the  input  axial  period  were  exact  then.  Figure  5.13  would  be  closed.  The 
stair  step  appearance  of  Figure  5.13  is  due  to  the  inputted  computer  step  size  --  the 
curve  is  actuallv  smooth 
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Figure  5. 1 3 

The  results  of  the  phasing  orbit  propagated  over  one  axial  period  for  a  6  X  6 
gravity  field,  with  drag,  are  shown  in  Figure  5.14.  Due  to  the  large  amount  of  data 
points  generated  for  the  6X6  gravity  field,  data  at  every  20th  ascending  nodal  passage 
was  plotted,  therefore,  the  appearance  of  the  graph  is  very  erratic.  If  data  were  not 
taken  at  every  20th  ascending  nodal  passage,  but  at  an  interval  consistent  with  Figure 
5  13.  then  Figure  5.14  would  appear  as  a  dark  mass  making  analysis  difficult.  The 
importance  of  Figure  5  14  is  not  in  its  erratic  shape,  but  like  Figure  5.13,  the  argument 
•  f  periapsis  and  the  eccentricity  both  are  bounded. 
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Figure  5.14 


The  MGCO  phasing  orbit  was  analyzed  to  gain  an  understanding  of  the  nature  of 
frozen  orbits.  As  defined  in  Chapter  1,  this  thesis  considers  a  frozen  orbit  as  any  orbit 
m  which  the  time  rate  of  change  of  one  or  more  of  the  orbital  elements  is  approximate^ 
zero,  or  nonsecular.  For  example,  the  above  orbit  (in  a  6  X  6  gravity  field)  does  not 
possess  a  single  orbital  element  whose  time  rate  of  change  is  zero;  however,  the 
argument  of  the  periapsis  oscillates  about  its  original  position,  and  hence,  the  phasing 
r b 1 1  is  considered  frozen.  Further  analysis  will  be  carried  out  to  determine  if  other 
frozen,  or  stable  orbits  exist  other  than  that  class  of  polar  orbits  with  the  periapsis 
located  over  the  poles.  Further,  are  there  orbits  whose  time  rate  of  change  of  one  or 
more  orbital  elements  equals  zero'1  If  so,  where  are  these  orbits  and  what  are  their 
ad  v  antages 1 


Semi  Ma n>r  Axis  Fqual  (2  4393.4  K ilometers 

Imtiallv.  a  value  of  the  semi  major  axis  of  4393  4  km  3nd  an  eccentricity  of  .1  was 
chosen  I  hose  values  establish  a  periapsis  altitude  of  approximately  560  km.  The  first 
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goal  is  to  freeze  one  of  the  orbital  elements.  -  ,  or  when  in  a  t>  X  6  gras  its  field 
1  r  the  MGCO  phasing  orbit,  the  argument  of  the  periapsis,  and  the  eccentricits  showed 
the  greatest  rate  of  change  oser  one  orbital  period,  therefore,  these  two  elements  will  be 
the  focus  of  this  step 

With  the  \alues  of  the  semi  major  axis  and  the  eccentricity  established,  the 
program  Omega  (see  Chapter  IV  and  Appendix  E)  was  run  in  order  to  determine  the 
salue  of  the  critical  inclination  angle  [that  \alue  of  inclination  that  "freezes"  ■<  oser  one 
orbital  period]  for  the  case  of  a  6  X  0  gras  its  field.  W  ith  this  salue  of  as  a  baseline,  a 
0  \  6  gras  its  field  was  introduced  and  numerous  runs  of  ASAP  were  made  in  order  t” 
find  a  critical  inclination  salue  of  68  15285662  degrees.  This  critical  inclination,  along 
with  the  other  orbital  elements  inputted  into  ASAP  define  Reference  Orbit  »l 
Reference  Orbit  *l’s  input  is  listed  in  Table  5.2 

Orbital  Elements  for  Reference  Orbit  **l 
Table  5.2 


Input 

Orbital 

Elements 

for: 

i  km 

** 

•  degrees 

i.  de¬ 
grees 

-  degrees 

v  de¬ 
grees 

Ref 

Orbit  #1 

4393  4 

1 

68.15285662 

90.00 

270  00 

90  00 

figures  5  15  through  5  17  indicate  that  the  change  in  <  for  this  orbit  is  indeed 
zero,  while  the  change  in  and  is  not  equal  to  zero  Further.  Figure  5  16  indicates 

that  the  change  in  the  semi  major  axis  vs  the  change  in  the  argument  of  periapsis  is 

bounded  Figure  5  18  shows  Reference  Orbit  # I  propagated  over  a  255  das  period 
Although  255  days  is  only  a  fraction  of  this  orbit's  axial  period  (the  axial  period  is  on 
the  order  of  15  sears),  it  is  sufficient  to  see  that  the  effect  of  the  change  in  eccentricity 

and  inclination  causes  the  argument  of  the  periapsis  to  change  by  approximately  200 

degree1-  This  does  not  compare  fasorably  with  the  MGCO  phasing  orbit 


58 


269.20 


269.70  270.20 

ARC  of  the  PERIAPSIS  (DEG) 


270  71) 


Arc  of  the  Periapsis  vs.  Inclination.  Ref.  Orbit  #1,  One  Orbital  Period.  6  ) 

Field 

Figure  5. 17 


turn  00  t  - 

■  i 


I)  111) 


ii.OO 


100.00 


PAYS 


200.00 


v  6  Gras  its 


000.00 


1  bunco  in  Arc  "t  the  Periapsis  Oser  255  Dass.  Ref  Orbit  »  1 .  6  \  6  Gras  its  1  icld 

f  igure  5  1 8 


I  he  above  figures  reveal  that  the  unbounded  nature  of  eccentricity  and  inclination 
adversely  effect  the  change  in  the  argument  of  the  periapsis.  The  analysis  indicates  the 
periupsis  does  not  oscillate  about  a  particular  point  (as  in  the  case  of  the  MGCO  phasing 
■  rbiti.  but  instead,  is  unbounded.  A  search  for  a  input  value  of  the  eccentricity  which 
causes  the  change  of  the  eccentricitv  and  the  inclination  to  be  zero  over  one  orbital 
period  was  made  for  values  of  eccentricitv  from  .01  to  0.7.  Figures  5.19  and  5.20  show 
the  results  of  this  search  and  reveals,  for  Reference  Orbit  »l,  a  value  of  eccentricitv 
which  drives  the  change  in  eccentricity  and  inclination  (over  one  orbital  period)  to  zero 
does  not  exist  Note,  any  eccentricity  greater  than  approximately  .227  will  cause  impact 
with  the  planet's  surface. 

F  igure  5  21  reflects  the  effects  of  various  eccentricities  and  inclinations  upon  the 
change  in  eccentricity  Since  circular  orbits  facilitate  the  use  of  scientific  instruments 
designed  to  observe  the  surface  of  \l3rs.  it  is  desirable  to  keep  the  value  of  the 
eccentricity  to  a  minimum  Also,  in  order  to  minimize  the  effects  of  atmospheric  drag  a 
minimum  periapsis  altitude  of  200  km  is  imposed.  Given  Figure  5.21  and  the  above 
restrictions,  analysis  revealed  that  an  eccentricity  of  0.3  and  a  semi  major  axis  of 
5133.428571  km  offers  the  best  compromise  between  the  desire  to  keep  eccentricity  to  a 
minimum,  and  the  need  for  an  eccentricity  which  drives  the  change  in  eccentricity  over 
me  orbital  period  to  zero 
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With  the  value  of  the  semi  major  axis  established  at  5133.428571  km,  the  value  of 
the  eccentricity  was  swept  from  e  =  0.01  to  <>  =  0.3  for  values  of  inclination  ranging  from 


1  to  90  degrees  (see  Figures  5.22  and  5.23). 


Inclination  vs.  Change  in  Eccentricity,  One  Orbital  Period,  Ref.  Orbit  #2,  6X6  Gravity 

Field 

Figure  5.22 

The  above  graph  shows  the  effects  of  various  values  of  eccentricity  and  inclination 
on  the  change  in  eccentricity  over  one  orbital  period.  From  this  graph  was  determined 
an  inclination  angle  of  15.05252881  degrees  that  will  cause  the  change  in  eccentricity  to 
equal  zero  over  one  orbital  period.  These  orbital  parameters,  along  with  the  other 
associated  input  parameters  define  Reference  Orbit  #2,  and  are  listed  in  Table  5.3. 


■T'  •  w\i  ~  V  m  \  -  V  j  *  U  9  ■'  *  V  *  T ■  V *  I  ■  'J  ■!1*';»'.|,V11 


In  Figure  5.23  the  effect  of  eccentricity  on  the  change  in  the  argument  of  the 
periapsis  over  one  orbital  period  is  investigated.  Three  important  findings  stand  out. 
First,  there  appears  to  be  values  of  eccentricity  near  zero  such  that  no  matter  what  the 
angle  of  inclination,  the  change  in  the  argument  of  the  periapsis  over  one  orbital  period 
will  never  equal  zero.  Second,  there  exist  values  of  eccentricity  and  inclination  (from  e 
=  0.03586336  at  i  =  90  degrees  to  <?  =  0.3  at  •  =  65.91286827  degrees)  which  cause  the 
change  in  the  argument  of  the  periapsis  to  equal  zero  over  one  orbital  period.  Third,  as 
the  eccentricity  increases  (at  least  from  0.03586336  to  0.3)  the  resulting  critical 
inclination  angle  decreases. 

The  value  of  the  angle  of  inclination  which  causes  the  change  in  the  argument  of 
the  periapsis  to  equal  zero  over  one  orbital  period  when  eccentricity  is  equal  to  0.3, 
together  with  the  other  orbital  inputs,  defines  Reference  Orbit  #3.  The  input  values  for 
Reference  Orbit  #3  are  listed  in  Table  5.3 


Orbital  Elements  for  Reference  Orbits  #2  and  #3 
Table  5.3 


Input 

Orbital 

Elements 

for: 

a  km 

e 

i  degrees 

n  de¬ 
grees 

w  degrees 

v  de¬ 
grees 

Ref 

Orbit  #2 

5133.428571 

.3 

15.05252881 

90.00 

270.00 

90.00 

Ref 

Orbit  #3 

5133.428571 

.3 

65.91286827 

90.00 

270.00 

90.00 

65 
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Table  5.4  list  the  input  orbital  elements  that  causes  the  change  in  the  argument  of 
the  periapsis  to  equal  zero  over  one  orbital  period  when  the  inclination  equals  90 
degrees.  This  orbit  is  known  as  Reference  Orbit  #4. 


Orbital  Elements  for  Reference  Orbit  #4 
Table  5.4 


Input 

Orbital 

Elements 

for: 

a  km 

e 

<  degrees 

o  de¬ 
grees 

w  degrees 

v  de¬ 
grees 

Ref 

Orbit  #4 

5133.428571 

.03586336 

90.00 

90.00 

270.00 

90.00 

Figures  5.24  through  5.26  show  Reference  Orbit  #2  over  one  orbital  period.  From 
these  three  graphs  it  can  be  seen  that  only  the  change  in  eccentricity  over  one  orbital 
period  is  zero.  Propagating  Reference  Orbit  #2  for  90  days  reveals  that  the  argument  of 
the  periapsis  changes  by  720  degrees  during  this  time  period  (see  Figure  5.27). 
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Figure  5.27 

Figures  5.28  through  5.30  show  Reference  Orbit  #3  over  one  orbital  period. 
Figure  5.31  shows  Reference  Orbit  #3  propagated  over  90  days.  Here  the  change  in  the 
argument  of  the  periapsis  over  a  90  day  period  is  only  34  degrees. 

Comparing  Reference  Orbits  #2  and  #3  shows  that  the  change  in  the  argument  of 
the  periapsis  and  the  semi  major  axis  are  significantly  larger  for  Reference  Orbit  #2. 
The  change  in  the  argument  of  periapsis  for  Reference  Orbit  #2  is  approximately  1.5 
degrees,  while  the  change  for  Reference  Orbit  #3  is  zero.  Likewise,  the  change  in  the 
semi  major  axis  for  Reference  Orbit  #2  is  approximately  0.5  kilometers,  compared  to 
approximately  zero  change  for  Reference  Orbit  #3.  The  situation  reverses  when  looking 
at  the  eccentricity  and  the  inclination.  The  change  in  eccentricity  over  one  orbital 
period  for  Reference  Orbit  #2  is  equal  to  zero,  while  the  change  in  eccentricity  for 
Reference  Orbit  #3  is  approximately  0.00005.  For  inclination,  Reference  Orbit  #3 
experiences  a  change  that  is  approximately  10  times  greater  than  that  of  Reference  Orbit 


Figure*  5.32  through  5.34  show  the  changes  for  one  orbital  period  associated  with 
Reference  Orbit  »4.  The  magnitude  of  the  change  of  the  argument  of  the  periapsis.  the 
eccentricity .  the  semi  major  axis,  and  the  inclination  are  all  of  the  same  order  as  those 
changes  for  Reference  Orbit  «3;  however,  for  Reference  Orbit  «4.  the  the  argument  of 
the  periapsis  changes  b>  150  degrees  per  90  days  (see  Figure  5.35),  as  opposed  to  34 
degrees  per  90  days  for  Reference  Orbit  #3.  The  reason  can  be  seen  in  Figure  5.22. 
For  Reference  Orbit  *3  the  inclination  is  increasing  with  each  orbit  causing  an 
increasing  smaller  change  in  the  eccentricity  for  each  successive  orbit.  For  Reference 
Orbit  «4  the  inclination  is  effectively  decreasing  with  each  orbit  causing  an  increasing 
larger  change  in  the  eccentricity  for  each  successive  orbit.  Figure  5.23  shows  that  an 
increase  in  eccentricity  and  inclination  (as  is  the  case  for  Reference  Orbit  #3),  and  an 
increase  in  eccentricity  associated  with  a  decrease  in  inclination  (Reference  Orbit  #4) 
both  induce  a  positive  change  in  the  argument  of  the  periapsis  over  one  orbital  period. 
Since  the  change  in  the  argument  of  the  periapsis  is  calculate  by  subtracting  the  final 
value  from  the  initial  value,  a  positive  change  implies  that  the  starting  value  for  the 
argument  of  the  periapsis  is  greater  than  the  value  of  the  argument  of  the  periapsis  one 
orbit  later;  therefore,  both  Reference  Orbits  #3  and  #4  are  experiencing  a  decrease  in 
the  argument  of  the  periapsis.  The  difference  in  the  magnitude  of  these  decreases  is  due 
to  the  initial  value  of  the  inclination  angle. 

For  Reference  Orbit  #3,  the  inclination  value  starts  out  being  the  critical 
inclination.  The  effect  of  the  increase  of  eccentricity  is  to  decrease  the  value  of  the 
critical  inclination  (see  Figure  5.23).  At  the  start  of  each  orbital  period.  Reference  Orbit 
»3's  inclination  has  increased  over  the  starting  value  of  the  previous  orbital  period  (see 
Figure  5.30).  The  combined  effect  is  that  Reference  Orbit  **3's  inclination  becomes 
slightly  greater  than  the  critical  inclination  angle.  This  causes  the  change  in  the 
argument  of  periapsis  over  one  orbital  period  to  be  slightly  positive,  thus  causing  the 
argument  of  the  periapsis  to  slowly  decrease 

Because  ot  the  initial  value  of  Reference  Orbit  #4’s  inclination  (  =  90  degrees),  the 
argument  of  the  periapsis  wants  to  decrease  at  its  maximum  rate  I  he  only  parameter 
holding  it  back  is  the  initial  eccentricity  This  eccentricity  increases  with  time,  and  with 


this  increase  in  eccentricity,  a  positive  change  (as  discussed  above)  in  the  argument  ot 
the  periapsis  occurs.  This  positive  change  in  the  argument  of  the  periapsis  enhances  the 
natural  tendenc>  for  an  orbit  of  this  inclination  to  decrease  the  argument  of  the 
periapsis  in  value.  Thus,  causing  the  change  in  the  argument  of  the  periapsis  to  be 
much  greater  than  that  of  Reference  Orbit  #3. 

At  this  value  of  the  semi  major  axis  and  eccentricity,  either  the  change  in  the 
argument  of  the  periapsis  or  the  change  in  the  eccentricity  over  one  orbital  period  can 
be  set  to  zero,  but  not  both.  Figures  5.24  and  5.27  indicate  that  selecting  an  inclination 
which  drives  the  change  in  eccentricity  to  zero  will  result  in  a  rapid  change  in  the 
argument  of  the  periapsis.  Thus,  in  the  effort  to  control  the  argument  of  the  periapsis, 
there  is  no  advantage  to  driving  only  the  change  over  one  orbital  period  in  the 
eccentricity  to  zero. 


vs  »,  Ref  Orbit  #4,  One  Orbital  Period.  6.X6  Gravity  Field 
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Figure  5.35 


Figure  5.23  shows  that  the  value  of  the  critical  angle  of  inclination  is  dependent 
upon  the  value  of  the  eccentricity .  However,  because  the  change  over  one  orbital  period 
of  the  eccentricity  is  non  zero  for  all  the  possible  values  of  the  critical  inclination  (see 
Figure  5.22),  the  value  of  the  critical  inclination  is  itself  changing  over  time,  thus 
inducing  a  change  in  the  argument  of  the  periapsis.  The  change  in  the  argument  of  the 
periapsis  is  the  slowest  at  the  maximum  allowable  eccentricity  U  =  0.3),  and  the  fastest  at 
the  minimum  allowable  eccentricity  <«  =  0.03586336)  .  The  question  arises,  for  *  =  0.3  is 
there  a  semi  major  axis  value  such  that  the  change  in  the  argument  of  the  periapsis,  the 
eccentricity,  the  semi  major  axis,  and  the  inclination  are  all  equal  to  zero?  If  so.  what 
are  the  characteristics  of  this  orbit,  and  what  effect  does  a  change  in  eccentricity  have 
upon  such  an  orbit0  The  next  part  of  this  analysis  will  focus  upon  these  questions. 


•Xnalvsis  From  5133  KM  Out  To  Geosynchronous 

Sweeping  the  value  ot  the  semi  major  axis  from  5133  km  to 
inclination  values  trom  I  to  90  degrees,  and  noting  the  change  over  one 


20.000  km,  for 
orbital  period  of 


75 


the  argument  ot  the  periapsis,  eccentricity,  inclination,  and  semi  major  axis  yields  the 
figures  shown  in  Appendix  F  through  H.  (Note,  for  Mars,  geosynchronous  occurs  at 
;n.400  km.)  Although  the  inclination  was  advanced  in  increments  of  5  degrees,  an 
increment  of  10  degrees  is  sufficient  to  show  the  trend,  and  is  used  in  presenting  these 
figures.  The  most  interesting  results  occur  at  an  inclination  of  approximately  70  degrees. 
Figures  5.. >6  through  5.38  highlight  these  results. 

At  a  semi  major  axis  value  of  approximately  17,000  km,  and  an  inclination  of 
approximately  70  degrees.  Figure  5.36  shows  the  change  in  the  argument  of  the  periapsis 
and  the  change  in  the  eccentricity  are  both  approximately  zero. 

Further  analysis  showed  that  the  zero  change  over  one  orbital  period  of  these  two 
parameters  (argument  of  the  periapsis,  and  eccentricity)  actually  occurs  when  the  semi 
major  axis  is  17,190  km  and  the  inclination  is  69.9750  degrees.  These  parameters, 
together  with  the  other  associated  input  parameters  define  Reference  Orbit  #5,  and  are 
shown  in  Table  5.5 

Orbital  Elements  for  Reference  Orbit  *5 
Table  5.5 
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The  time  rate  of  change  of  the  argument  of  the  periapsis,  the  eccentricity,  the  semi 
major  axis,  and  the  inclination  for  Reference  Orbit  #5  are  shown  by  Figures  5.39 
through  5.41.  Figure  5.42  reveals  that  although  the  change  in  the  argument  of  the 
periapsis  is  zero  over  one  orbital  period,  over  an  extended  time  the  argument  of  the 
periapsis  decreases  at  the  rate  of  1.15  degrees  per  90  days.  This  is  because  of  the 
combined  effect  of  the  change  in  the  semi  major  axis,  on  the  order  of  0.1  km  per  orbit 
(see  Figure  5.40),  and  the  0.00026  degree  per  orbit  change  in  the  inclination  (see  Figure 
5.41).  From  Figure  5.36  it  can  be  seen  that  a  small  decrease  from  the  semi  major  axis 
value  of  17,000  km  will  induce  a  slow  decrease  in  the  eccentricity  and  the  argument  of 
the  periapsis.  Comparing  Figures  5.36  and  F.9  (see  Appendix  F)  reveals  an  increase  in 
inclination  will,  at  this  particular  value  of  the  semi  major  axis,  also  result  in  a  slight 
decrease  in  the  argument  of  the  periapsis  over  time.  Figures  5.36  through  5.38  show 
that  as  the  semi  major  axis  slowly  decreases,  the  change  in  the  argument  of  the 
periapsis,  the  change  in  the  eccentricity,  the  change  in  the  inclination,  and  the  change  in 
the  semi  major  axis  will  progressively  increase.  Due  to  Reference  Orbit  #5’s  long 
orbital  period  (approximately  20  hours)  a  0.1  km  decrease  in  the  semi  major  axis  per 
orbit  results  in  only  an  approximate  44  km  decrease  in  the  semi  major  axis  per  year. 
Hence,  although  the  change  in  the  above  orbital  elements  occurs  at  a  progressively 
increasing  rate,  the  increase  in  rate  is  painstakingly  slow.  This  accounts  for  the  constant 
slope  of  Figure  5.42. 

Equation  (5.1)  shows  the  predominant  effect  of  an  increase  in  the  semi  major  axis 
is  a  decrease  in  the  rate  of  change  of  the  argument  of  the  periapsis.  Therefore,  the  slow 
rate  of  Reference  Orbit  #5's  change  in  the  argument  of  the  periapsis,  when  compared  to 
Reference  Orbit  #3,  is  not  due  entirely  to  any  special  effects  of  one  change  in  an  orbital 
parameter  cancelling  the  effects  of  the  change  in  another  orbital  parameter.  Figure  5.43 
shows  the  change  in  the  argument  of  the  periapsis  for  an  orbit  of  the  same  semi  major 
axis  as  Reference  Orbit  #5,  but  at  45  degrees  inclination.  A  comparison  of  Figures  5.42 
and  5.43  reveals  Reference  Orbit  #5  experiences  the  same  magnitude  of  change  in  one 
>ear  as  the  change  experience  over  a  90  day  period  for  the  orbit  in  Figure  5.43.  Thus 
indicating  Reference  Orbit  #5  is  the  more  stable  orbit. 
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Figure  5.43 

Appendices  I  through  K  show  the  effects  of  varying  the  eccentricity  and  semi 
major  axis  upon  the  change  in  the  argument  of  the  periapsis,  the  eccentricity,  the 
inclination  angle,  and  the  semi  major  axis.  Throughout  these  figures  the  input 
inclination  angle  remains  70  degrees.  These  appendices  offer  trend  information,  and 
because  not  atl  combinations  of  eccentricity  and  semi  major  axis  exist  without  causing 
impact  with  the  planet  Mars,  caution  must  be  exercised  in  using  these  figures.  From 
Appendix  I  it  can  be  seen  that  for  eccentricities  from  0.01  to  0.3  (see  Figure  5.36),  the 
change  in  the  argument  of  the  periapsis  over  one  orbital  period  has  its  zero  value 
between  approximately  17,000  and  18,000  kilometers.  For  eccentricities  between  0.01 
and  0.6  the  change  in  eccentricity  will  also  become  zero  somewhere  between  17,000  and 
18,000  km.  As  previously  mentioned,  only  when  the  eccentricity  is  0.3  is  there  one 
value  of  the  semi  major  axis  that  simultaneously  drives  both  values  to  zero.  Further, 
between  the  semi  major  axis  values  of  approximately  12,000  and  13,000  km  there  is 
another  region  where  the  change  in  eccentricity  over  one  orbital  period  becomes  zero. 
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Appendix  J  shows  that  at  an  eccentricity  of  0.01  (Figure  J.l)  the  change  in 
inclination  is  fairly  insensitive  to  changes  in  the  semi  major  axis  from  approximately 
13,000  to  17,000  km.  As  eccentricity  increases.  Figures  J.l  through  J.6  show  that  the 
change  in  the  inclination  becomes  more  sensitive  to  the  value  of  the  semi  major  axis; 
however,  although  not  exactly  zero,  the  change  in  the  inclination  over  one  orbital  period 
remains  relatively  constant,  and  approximately  zero  in  the  semi  major  axis  region  of 
17,000  to  18,000  km.  Also,  starting  at  an  eccentricity  of  approximately  .1,  the  change  in 
inclination  over  one  orbital  period  is  approximately  zero  between  semi  major  axis  values 
of  approximately  12,500  and  13,000  km. 

In  Appendix  K,  the  change  in  the  semi  major  axis  over  one  orbital  period  appears 
to  be  fairly  sensitive  to  changes  in  both  the  semi  major  axis  and  the  eccentricity. 
Through  out  the  range  of  values  of  the  eccentricity  there  appears  two  regions  where  the 
zero  change  in  the  semi  major  axis  exist.  These  regions  exist  from  a  semi  major  axis  of 
approximately  12,200  to  13,200  km,  and  approximately  16,000  to  18,000  km. 

Atmospheric  Drag 

In  the  next  phase  of  this  analysis,  atmospheric  drag  was  introduced  to  the  above 
orbits.  In  all  cases  the  atmospheric  drag  showed  no  appreciable  effect  over  one  orbital 
period.  Reference  Orbit  #3  was  propagated  over  a  one  year  period,  both  with  and 
without  atmospheric  drag.  The  results  show  that  when  in  the  presence  of  drag,  the 
eccentricity  decreased  by  0.00002986,  and  the  semi  major  axis  decreased  by  798  meters 
more  than  if  atmospheric  drag  were  not  present.  Reference  Orbit  #5  was  also 
propagated  over  a  one  year  period.  Here  the  results  showed  no  appreciable  effects  when 
in  the  presence  of  drag.  This  finding  is  not  surprising  given  that  the  periapsis  altitude 
of  Reference  Orbit  #5  is  approximately  8,600  km.  Because  of  the  height  of  the  orbits 
investigated,  atmospheric  drag  effects  are  minimal. 


Conclusions 


From  the  analysis  section,  two  general  classifications  of  results  were  found.  The 
first  involves  the  characteristics  of  the  orbital  parameters  effecting  the  control  of  the 
argument  of  the  periapsis,  and  the  second  involves  the  two  regions  of  relative  orbital 
stability. 

Characteristics  of  the  Orbital  Parameters  Effecting  the  Control  of  thg.  Argument  of 
the  Periapsis.  For  the  MGCO  phasing  orbit  none  of  the  orbital  elements  (u>,  <>,  a,  and  t) 
experienced  a  zero  time  rate  of  change  over  one  orbital  period  when  in  the  presence  of  a 
6X6  gravity  field.  Yet,  the  values  of  the  argument  of  the  periapsis  and  the 
eccentricity  are  bounded.  The  first  step  of  the  analysis  sought  to  discover  the  nature  of 
the  argument  of  the  periapsis  when  the  time  rate  of  change  of  the  argument  of  the 
periapsis  is  zero  over  one  orbital  period.  Through  a  careful  selection  of  the  inclination 
angle,  the  change  in  the  argument  of  the  periapsis  was  driven  to  zero  over  one  orbital 
period.  The  results  showed  that  a  zero  change  in  the  argument  of  the  periapsis  has 
associated  with  it  a  non  zero  change  in  the  eccentricity  and  the  inclination  angle  (see 
figures  5.15  through  5.18).  These  two  non  zero  parameters  induce  a  long  period  change 
in  the  argument  of  the  periapsis  that  is  not  bounded,  but  rather  periodic.  Further,  from 
Figure  5.23  it  can  be  seen  that  the  critical  inclination  angle  has  a  range  of  values, 
dependent  upon  the  eccentricity  of  the  orbit.  From  Figures  5.31  and  5.35  it  is  revealed 
that  the  rate  of  change  in  the  argument  of  the  periapsis  that  is  induced  by  the  change  in 
eccentricity  and  inclination  is,  as  in  the  case  of  a  6  X  0  gravity  field,  very  sensitive  to 
the  initial  value  of  the  critical  angle  of  inclination.  If  the  eccentricity  is  such  that  a 
critical  angle  of  inclination  has  a  value  that  is  close  to  90  degrees,  the  rate  of  change 
induced  in  the  argument  of  the  periapsis  will  be  much  greater  than  the  case  where  the 
critical  inclination  is  near  some  lower  value  of  inclination.  Thus,  driving  only  the 
change  in  the  argument  of  the  periapsis  to  zero  is  not  sufficient,  when  in  the  presence 
of  a  6  X  6  gravity  field,  to  control  the  argument  of  the  periapsis. 


In  the  next  step  of  the  analysis,  a  value  of  the  semi  major  axis  and  the  inclination 
was  selected  that  allows  the  time  rate  of  change  over  one  orbital  period  of  the 
eccentricity  to  be  driven  to  zero.  Figure  5.24  shows  that  for  this  orbit  there  is  a  large 
change  over  one  orbital  period  in  the  argument  of  the  periapsis;  hence,  driving  the  time 
rate  of  change  in  the  eccentricity  to  zero  will  not  result  in  the  desired  bounded  condition 
of  the  argument  of  the  periapsis. 

Searching  values  of  the  semi  major  axis  ranging  from  5133  km  to  20,000  km,  for 
inclinations  ranging  from  1  to  90  degrees  lead  to  the  discovery  of  an  orbit  in  which  both 
the  argument  of  the  periapsis  and  the  eccentricity  are  zero  over  one  orbital  period. 
Figure  5.39  shows  that  the  argument  of  the  periapsis  and  the  eccentricity  are  indeed 
bounded;  however,  the  semi  major  axis  and  the  angle  of  inclination  are  not  bounded. 
Figures  5.36  through  5.38  indicate  how  the  unbounded  nature  of  the  semi  major  axis  and 
the  inclination  effect  the  argument  of  the  periapsis  and  the  eccentricity.  The  results 
indicate  that  in  the  presence  of  a  6  X  6  gravity  field,  control  of  the  argument  of  the 

periapsis  is  not  gained  by  driving  the  short  term  perturbations  in  the  argument  of  the 

periapsis  and  the  eccentricity  to  zero. 

Regions  of  Relative  Orbital  Stability.  From  this  analysis  it  is  evident  that  driving 
the  short  term  perturbations  of  one  or  two  of  the  orbital  elements  is  not  sufficient  to 
control  the  argument  of  the  periapsis.  Rather,  the  short  term  perturbations  for  the 
change  in  the  argument  of  the  periapsis,  the  eccentricity,  the  semi  major  axis,  and  the 
angle  of  inclination  must  all  be  driven  to  zero.  Appendices  F  through  K  show  that  at 
eccentricities  from  0.01  to  0.6  an  orbit  that  freeze  the  argument  of  the  periapsis,  the 
eccentricity,  the  semi  major  axis,  and  the  inclination  does  not  exist.  The  best  that  can 
be  obtained  is  that  the  change  in  three  out  of  four  of  the  orbital  elements  can  be  driven 

to  approximately  zero.  Appendices  F  through  K  indicate  that  this  takes  place  in  two 

distinct  regions.  The  first  being  for  a  semi  major  axis  from  approximately  17,000  km  to 
18,000  km,  with  an  eccentricity  ranging  from  approximately  0.01  to  0.6  and  an 
inclination  value  of  approximately  70  degrees.  Within  this  region  the  change  in  the 
argument  of  the  periapsis,  the  change  in  the  eccentricity,  and  the  change  in  the 
inclination  can  be  driven  to  approximately  zero,  while  the  change  in  the  semi  major  axis 


prominently  remains  non  zero.  The  second  region  exist  for  semi  major  axis  values  from 
approximately  12,000  km  to  13,000  km,  with  an  eccentricity  ranging  from  0.01  to  0.6 
and  an  inclination  value  of  approximately  70  degrees.  Within  this  region  the  change  in 
the  eccentricity,  the  semi  major  axis  and  the  inclination  can  be  made  approximately 
zero,  but  the  change  in  the  argument  of  the  periapsis  can  not  be  driven  to  approximately 
zero. 

Atmospheric  Drag..  At  the  beginning  of  this  research  it  was  thought  that  the 
locations  for  the  frozen  and  stable  orbits  that  might  be  found  would  be  near  the  planet’s 
surface.  This  was  not  the  case.  In  fact  the  orbits  looked  at  were  of  sufficient  height 
that  the  atmospheric  drag,  even  when  propagated  over  a  one  year  period  had  only  very 
slight  effects  on  the  semi  major  axis,  and  no  discernible  effects  on  the  eccentricity. 

Recommendations 

Examining  the  MGCO  phasing  orb«t  reveals  that  the  change  in  the  argument  of  the 
periapsis  and  the  inclination  are  both  negative  (values  increase  over  one  orbital  period), 
while  the  change  in  the  eccentricity  and  the  semi  major  axis  are  both  positive  (values 
decrease  over  one  orbital  period).  Appendices  F  through  H  show  that  for  an  eccentricity 
of  0.3,  a  region  exist  from  a  semi  major  axis  of  approximately  13,000  km  to  17,100  km, 
and  an  inclination  value  of  approximately  35  degrees  to  65  degrees  where  the  same 
characteristics  of  change  in  the  argument  of  the  periapsis,  eccentricity,  semi  major  axis, 
and  the  inclination  exist  as  exist  for  the  MGCO  phasing  orbit.  The  next  step  in  any 
follow  on  study  ought  to  focus  on  this  region.  If  this  region  does  prove  to  have 
bounded  changes  in  the  argument  of  the  periapsis,  then  further  analysis  needs  to  be 
made  at  other  values  of  the  eccentricity. 

It  is  interesting  to  note  that  the  above  region,  and  the  regions  of  relative  orbital 
stability  described  earlier  occur  between  the  moons  of  Mars,  Phobos  (mean  distance  of 
9,380  km)  and  Deimos  (mean  distance  of  23,474  km).  Although  the  mass  of  these  moons 
are  slight,  they  will  have  a  perturbative  effect  upon  the  regions  of  relative  orbital 


Let  the  entire  mass  of  the  planet  exist  as  a  point  in  space,  then,  surround  this 
point  with  a  "simple"  surface  S,  The  surface  S  is  called  simple  if  it  has  a  finite  area  and 
does  not  have  points  that  intersect  or  touch  other  points  on  this  surface,  (see  Figure  Al) 
Each  small  area  of  the  surface  (da)  will  have  an  associated  normal  vector  n.  The 
procedure  is  to  determine  how  much  of  the  acceleration  (due  to  the  mass  v)  is  along  the 
normal  of  each  infinitesimal  area  of  surface,  and  then  to  integrate  over  the  entire  area 
of  the  surface  S.  This  will  yield  the  amount  of  acceleration  which  mass  v  exerts  over 
the  entire  surface  S  (19:49). 


Mathematically  this  process  is  modeled  as: 


fa  ■  n  da  =  <fc  -  r  •  ft  da  =  -CM  4  f-  •  n  da 


-  s  r 


(  a  i .  i : 


r  r 


Since  S  is  a  simple  surface,  the  value  of  the  above  integral  will  not  be  dependent 
upon  the  size  of  the  surface.  Therefore,  let  S  be  a  unit  sphere.  Then 


p  a  ■  nda  =  - CM  f  da  =  -CM  4 n r 

'  s  *  j  s 


=  -  CM  4  n 


Employing  the  Divergence  Theorem  of  Gauss  (12:440): 


(A  1 .2) 


a g ■  nda 


-f,7°d  *' 


equation  (A  1.2)  becomes: 


j>  a  g-  nda  =  J*  7  •  a  Qd\'  =  -  4nC M 


(A  1 .3) 


(A  1 .4) 


where  i  equals  the  volume  enclosed  by  the  surface  S. 

The  above  equation  assumes  the  entire  mass  of  the  planet  exists  as  a  point  mass 
located  at  the  center  of  a  unit  sphere.  This  restriction  is  removed  by  assuming  that  the 
mass  of  the  planet  is  evenly  distributed  throughout  the  unit  sphere  by  allowing: 


-I 


m  =  pd  1' 


where  o  equals  the  density  of  the  mass.  Equation  (A  1.4)  becomes: 


(A  1 .5) 


j  ■  afdl' = -4nG  j  pdl' 


(A  1 .6) 


90 


Because  equation  (A  1.6)  is  independent  of  the  size  or  shape  of  the  volume,  equate  the 
integrands  to  obtain: 


but 


V  •  a 5  =  -4 nCp 


(A  1 .7) 


so  equation  (A  1.7)  becomes: 


a,  =  ?  V 


(A  1 .8) 


7-  VI/  =  -4/rCp  (A  1 .9) 

V2F  =  -4 nCp 

Equation  (A  1.9)  is  known  as  Poisson’s  Equation.  This  equation  is  only  valid  for 
regions  within  the  planet’s  interior  (5:108).  Since  the  satellite  will  be  orbiting  outside 
the  planet’s  surface,  p  becomes  zero,  and  equation  (A  1.9)  becomes: 

V2l/  =  0  (A  1 .10) 

Equation  (A  1.10)  is  Laplace’s  Equation. 


Identities 


cos,  a  +  b  =  cos  a  cos  b  -  si n  a  sin  b 

sin  a+b!  =  sinacosb  +  cosasinb 

cosia-b  —  cosacosb  +  sinasinb 

sin  a-bl  =  sinacosb-cosasinb 


1  , 


cosacosb=-  cosici  +  b  +  cos  a-b)l 
2 


1  • 


sin  a  sin  b  =  ^  cos'  a  -  b  -  cos!  a  +  t> ) ] 


1  , 


s/na  cos  b  =  -  i  sin  '  a  +  b  ♦sinia-6-' 

2' 


cos  a  sin  b=-sin!a  +  b)-sin(a-6)] 


Euler’s  equations  are: 

eia -e  ia 

sin  a  =  — - - 

2/ 

eia  +  e 

cosa  =  — — - 


where  )  = 

e'a  =  cosa  +  /sin  a 


Binomial  Expansions  q£  £21  Qix  and  sin  mx 


Let  (3:2) 


rosmv*  real  part  <  e,m'  -  RE  e‘ 


(B2. 1 ) 


cos  mx  =  RE  e' 


cosmx ■=  RE  icos.v  +  /sin.v/' 


Noting  that  (1:11) 


-  b  n  =  > 


(B2.2) 


(B2.3) 


(B2.4) 


Equation  (B2.3)  can  be  written  as: 


■ RE  £(?)'■• 


(B2.5) 


sin  m  v  =  RE  [  -  ye' 


Then,  equation  (B2.6)  becomes: 


(B2.6) 


sinrn.v=  RE  (  y 1  1  cos"  ’  x  sin  *  .v 
« .  ft  V  s  J 


sions  of  Sin  »  mx  Cos  b  mx 


(B2.7) 


Multiplying  equtions  (B1.9),  and  (B1.10)  yields: 


„  „  [  e'^-e  7  " 

sin  mx  cos  mx  =  -  - — - 

V  2  J  J  \  2 


(B3.I) 


sin  mx  c  os  mx 


2 


p 


i 


> 


x 


b 


jxd 


which  becomes 


Q  b 

sin  .v cos  v  = 


a 

c 


b 

d 


c  q  /x  l  a  *  b  -  2  c  -  2  d  ' 


where 


p''*1'0  2c  2d  =  cos  a  *  b  -  2c  -  2d  v»ysinn  +  b-  2c-2d'.v 
In  equation  (B3.3),  let  a-t-m-2<-s  and  t>-m-s  (Born:5).  Then  (B3.3)  becomes: 


/ x  ci  *  b  2 f  2d 

P 


=  cos  L  -  2t  -  2c  -  2d  x  *  j  sin'  L  -  2t  -  2c  -  2d  .v 


(B3.4) 


Table  C.l  cont.  The  Inclination  Function 


Table  C.l  cont.  The  Inclination  Function 


Table  C.2  The  Eccentricity  Funtion  ( 10:38) 


Let  the  mass  of  the  unit  volume  of  atmosphere  be  denoted  by  m.  Then  (7:2): 


I-v.-v, 

m  =  - = - 

x>. 

where  m  =  the  mean  molecular  mass  of  the  unit  volume  atmosphere 

v.  =  the  number  of  molecules  of  type  i  per  unit  volume 

w.  =  the  molecular  mass  of  type  i  molecule 


Since  the  unit  volume  element  of  the  atmosphere  is  stationary  (i.e.  no  thermal  effects  are 
being  considered): 


f  IP  +  ^  DOWN  =  0 


(D1.2) 


where 


Fv  P  =  ALPlz+dz)-P  (=!] 
F  now s  =  y~  -S’  i  \f ,){  Adz  \g 


Applying  equations  (D1.3)  to  (D1.2)  yields: 


(D1.3) 


£ 


P  z*dz  -P  r.  1,)\dz]g  (D1-4) 

dP  =  -g  X'V.V.ic/s 

where  p  =  the  total  pressure  of  the  atmosphere  at  height  z 

i  =  the  area  on  the  surface  of  the  planet  subtended  by  the  unit  volume  of 

atmosphere 

/  =  the  acceleration  due  to  gravity 


From  the  ideal  gas  law  an  expression  for  the  total  pressure  can  be  written  as  (9:12): 

T \,;RT  (D1.5) 

D  _  • 


where  p  =  the  universal  gas  constant 

!  =  the  temperature 


=  the  volume 


mmmm 


=  p  =  density 


So  equation  (D1.5)  becomes: 


P  =  pRT 


Writing  equation  (D1.4)  in  terms  of  equation  (D1.7)  yields: 

p  *  cl z  '  RT  =  -g ,  £\v,  M,  )dz  (D1.8) 

p  z  +  dz,  - - dz 

H  RT 

dp - w~d~ 

Noting  that  for  the  unit  volume  element,  r-i,  dividing  equation  (D1.8)  by  equation 
(D1.6),  and  applying  equation  (Dl.l)  yields: 


dp  _  _ 

P  Y.N‘RT 


(D1.9) 


Taking  the  integral  of  both  sides  of  equation  (D1.9)  yields: 


.  grn 

n  p  =  - — = z  +  c 
H  RT 


(D1.10) 


When  -o,  r-inp..  As  a  result  equation  (DI.IO)  becomes: 


In  ^  =-|? 

\Po)  Rr 


(Di.i  i : 


Taking  the  exponential  of  both  sides  of  equation  (Dl.l  1 )  yields: 


I  Sffn  \ 

P-Poexp\-~^ 


(DM2) 
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C234567* 


This  program  calculates  the  geopotential  field  around 
the  planet  Mars.  It  writes  data  into  file  MARSGRAV.DAT 
in  three  columns  (corresponding  to  X,Y,Z  coordinates) 
where 

X  =  Longitude 
Y  =  Latitude 

Z  =  the  value  of  the  Geopotential 


C 
C 
C 
C 
C 
C 
C 
C 

CUritten  by  J.  W.  Foister,  III 
COate  10  Sept  87 

C 

C234567*********** ***************** ***************************** 

SNOFLOATCALLS 

PRfinRAM  MARS1 

REAL*8  C(19,19),S(19,19),LP(19, 19),LAT,PHI ,PI , 

^1  MU,RP,R,LONG,V,PT(91,91),B 

C  C  =  The  nondimensional  c  coefficient  to  the  grav  model  * 

C  S  *  The  nondimensional  s  coefficient  to  the  grav  model  * 

C  LP  =  The  values  of  the  Legendre,  and  associated  Legendre* 

C  polynomials  * 

C  LAT  s  The  latitude  * 

C  PH(  =  The  sin  of  the  latitude  * 

C  PI  *  The  classic  irrational  ncrnber  * 

C  MU  «  The  universal  gravitational  constant  multiplied  by  * 

C  the  mass  of  the  planet  Mars  * 

C  RP  =  the  equatorial  radius  of  Mars  * 

C  R  =  The  distance  from  the  center  of  Mars  at  which  it  is* 

C  desired  to  calculate  the  geopotential  * 

C  LONG  =  The  longitude  * 

C  V  =  An  intermediate  value  of  the  geopotential  (calculat* 

C  ion  not  yet  complete)  * 

C  PT  =  The  value  of  the  geopotential  at  a  particular  point* 

C  B  =  Intermediate  step  in  calculating  geopotential  * 

C234567********************************************************** 
C 
C 
C 

INTEGER  i , j ,k,m,n 

C234567*****  Determine/Set  the  intial  values  ******************** 
PI =4 .00*0 AT AN ( 1 .00) 

MU=4. 282828704 
RP=3393.400 
R=3893.400 
LAT=-90.D0 
LONG=0.D0 
V=0.00 
C234567 

OPEN ( 1 ,  F  ILE=  ‘MGRAV  ,  STATUS= 1  OLD '  ) 

C234567  This  file  contains  an  18  by  18  gravity  model  of  Mars 
C  Source:  Jet  Porpulsion  Laboratory,  EM  312/87-153 

C  20  April  1987 

C234567 

cci , i  )=i  .oo 

S(  1,1 )=0.00 
C(2, 1 )=0.00 
S(2, 1 )=0.00 
C<2,2)=0.00 
S{ 2 , 2 ) =0 .00 
00  10  i =3 ,19 
00  20  j  =  1 , i 

READ( 1 , ' (25X, F15. 13,15X,F15.13>’)  C( i , j ),S( i , j ) 

20  CONTINUE 

10  CONTINUE 
C234567 


Data  File  MGR  A  V.  Mars  Gravity  Model 

The  following  is  an  18  by  18  gravity  model  of  Mars  (see  reference  13). 
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.2800576750-14 
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.4087416430-16 
.7245096500-06 
.66544164  70-06 
.4128146200-07 
.2285234860-08 
.1389907480-09 
.4301416080  10 
.2814623440-11 
.4225581350-13 
.6068102800-14 
.4016611640-14 
.1844969620-15 
.4765389990-16 
.2957768750-17 
.5431980320-05 
.1055102150-08 
.1115885050-07 
.3135528130-08 
.1418926900-09 
.2135356910-10 
.7581287240-12 
.1937441730-12 
.8892082760-14 
.1618950470-14 
.2480011540-16 
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.4301486450-11 
.5252841020-12 
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.2103417510-05 
.2410731710-06 
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1  C234567********************************** *************************** 

2  C  Omega  is  a  program  that  finds  the  roots  to  delta  omega  * 

3  C  where  the  equations  for  delta  omega  (change  in  arg.  of  the  * 

4  C  periapsis)  are  found  in  THE  MOTION  OF  A  SATELLITE  IN  AN  AX  I  -  * 

5  C  SYMMETRIC  GRAVITATIONAL  FIELD,  by  R.  H.  MERSON.  As  found  in  * 

6  C  the  Geophysical  Journal  Vol  4,  1961,  p.17.  This  program  finds* 

7  C  the  roots  of  delta  omega  as  a  function  of  f,  f=(sin  i)**2.  * 

8  C  The  primary  equation  is:  delta  omega  =  cf**3*bf**2*af»d=0.  * 

9  C  Where  a,b,c,  and  d  are  constants  depending  on  the  semi  major  * 

10  C  axis  and  the  eccentricity  * 

11C  * 

12  C  WRITTEN  BY  J.W.  FOISTER,  III  * 

13  C  DATE  AUG  5  1987  * 

14  C  * 

15  C  J(1)  =  * 

16  C  J(2)  =  C20  coefficient  * 

17  C  J(3)  =  * 

18  C  J(4)  =  C40  coefficient  * 

19  C  J(5)  =  * 

20  C  J(6)  =  C60  coefficient  * 

21  C234567************************************************************* 

22  SNOFLOATCALLS 

23  PROGRAM  OMEGA 

24  C 

25  C 

26  C234567**********  DEFINE  THE  VARIABLES  **************************** 

27  REAL'S  J(6),LAT,A,B,C,D,E,F(3),EC,DLTA,P,0,R,X(3),Y,THETA(3), 

28  1  SEMI ,PI ,W(3),TP,AHLD 

29  C 

30  C 

31  C234567 

32  INTEGER  i.k, COUNTER 

33  C 

34  Pl=4.00*DATAN(1.00) 

35  R=1 .DO 

36  C234567*******  fo( lowing  values  are  for  the  planet  Mars  ******** 

37  J(2)=-0. 1960454460-02 

38  J(4)=0. 1889436780-04 

39  J(6)=0. 1340756980-05 

40  C 

41  OPEN ( 1 , F I LE= 1 LPT1 1 ) 

42  200  FORMAT ( IX, 1 f  =  « , F30. 15) 

43  210  FORMAT (1X,'f',i1,'  =  ■ , F30. 15 ) 

44  C 

45  C234567*******  INPUT  THE  DATA  *********************************** 

46  WR I TE(* , 100) 

47  100  FORMAT(1X, 'Semi  Major  axis  equals _ in  KMs. . .(F30. 15) ....', \) 

48  READ<*, 1 (F30.15) ' )  SEMI 

49  WR I TE( 1 , 105 )  SEMI 

50  105  FORMATf IX, ' The  Semi  Major  axis  in  Km  is  1 , F30 .15) 

51  WR I TE(*, 105)  SEMI 

52  C234567  convert  from  KM  to  Du's  this  is  for  a  Mars  orbit 

53  SEMI=SEMI/3393.4 

54  C234567 

55  WRI TE(*, 1 10) 

56  110  FORMAT(1X, 'Eccentrici ty  equals _ (F30.15) . 1 , \ ) 

57  REA0(* , ' (F30. 15) ' )  EC 

58  AHLD=(3.DO/2.DO) 

59  TP=(2*PI)*(SEMI**AHL0) 

60  C234567******  This  converts  from  Mars  Time  Units  to  Minutes  ******* 

61  TP=TP*( 15. 9197403800) 

62  C 

63  WR1TE(*,120)  SEMI, TP, EC 

64  WRITE(1,120>  SEMI, TP, EC 

65  120  F0RMAT(1X,  'For  the  Semi  Major  axis  *  \F30.15, 

66  1  '  the  period  (in  minutes)  is  '.F30.15,/, 

67  2  '  and  Eccentricity  *  '.F30.15) 

68  C234567*********  CALCULATE  THE  VALUES  OF  THE  COEFFICIENTS  ********** 

69  LAT =SEMI *( 1 -EC**2> 


I, 


C  SEMI  =  semi  major  axis,  lAT  -  semi  (atus  rectun 
J(1)=J(4)*((R/LAT)**4) 

J(3)=J(6)*((R/LAT)**6) 

J(5)=(J(2)**2)*((R/LAT)**4) 

C  A  through  0  are  coefficients  described  in  line  8  of  this  program 
A  =  (  (15.D0/4.D0)*J(2)*((R/LAT)**2)) 

A=A*J ( 1  )*( (930. DO/32. DO )*(945. DO/32. D0)*EC**2) 

A=A*J(3)*(( -33600. D0/320. DO)  (( 22575 . DO/64 .DO )*£C**2 ) 

1  (< 14175. DO/128. D0)*EC**4>) 

A=A*J(5)*((8S5.D0/48.D0)  ( (27 . DO/32 . D0)*EC**2 ) > 

C234567 

B=J(1 )*(( -735.00/32. DO)  ((2835. DO/128 ,D0)*EC**2 ) ) 

B=B*J( 3)*( (541800. 00/2560. 00 )♦< (343350. D0/512. D0)*EC**2) 

1  *((51975. 00/256. 00)*EC**4)) 

B=B*J(5)*((  4005.00/192.00) •(( 135 .DO/128. D0)*EC**2) ) 

C234567 

C=J(3)*(( ■ 1247400.00/10240.00)  ((381 150 .D0/5 1 2 .D0)*EC**2 > 

1  ((225225. DC/2048. 00)*EC**4)) 

C234567 

D=3.D0*J(2)*((R/LAT)**2)*J(1>*(( -240.00/32.00) 

1  ( (270. 00/32. D0)*EC**2))*J(3)*( (4200. D0/320. DO) 

2  ♦(( 3150.00/64.00 )*EC**2)*(( 1050. DO/64 .DO )»EC**4 ) ) 

3  *J(5)*((63. 00/48. D0)*EC**2) 


93  C234567 


94  P=((A/C)((B/C)**2)/3.D0) 

95  Q=((2.D0/27.00)*((B/C)**3)  ( ( A*B )/(3 ,00*C**2) )*(D/C) ) 

96  C234567 

97  0LTA=( 27.00*0**2) (4.00*P**3) 

98  C 

99  If  (DLTA  .  LT .  0.0)  THEN 

100  C  one  real  root  exist 

101  GOTO  500 

102  ELSE  1 F  (OLTA  .GE.  0.0)  THEN 

103  C  all  roots  are  real,  if  dlta  =  0  then  two  of  the  roots  are  the  same 

104  c  if  dlta  >  0  then  all  three  roots  are  different 

105  GOTO  700 

106  ELSE 

107  END  I F 

108  C234567*******  NEUTON  -  RAPHSON  METHOO  »•***•****•»*«»*»«***♦***** 

109  500  COUNTERS 

110  C234567set  inital  guess  of  f  value 

111  X ( 1 ) =0 . 500 

112  503  IF  (COUNTER  .GT.  100)  THEN 

113  UR  1 TE( 1 , 505 )  COUNTER 

114  505  FORMAT ( IX , ‘After  1 , i 3 , 1  iterations') 

115  F(1  )=X(2> 

116  GOTO  800 

117  ELSE 

118  END  I F 

1 19  U(1  )=(2.D0*C*X(1 )**3)*(B*X(1>»*2)0 

120  U(2)=(3.00*C*X( 1 )**2)*(2.00*B*X( 1 )  )*A 

121  IF  (U(2)  .EQ.  0.0)  THEN 

122  UR  1  TEC  1 , 510) 

123  510  FORMAT ( lx ,' F i rst  derivative  =  0,  program  stopped1) 

124  GOTO  900 

125  ELSE 

126  END  IF 


127  C234567 


X(2)=W(1)/U(2) 
x(3)=x(2)  X(1) 

X(3)=DABS(X( 3) ) 

IF  (X(3)  .LE.  1.D-12)  THEN 
C the  root  to  the  equation  is  X(2) 

F ( 1 ) =X( 2 ) 

COUNTERS 
GOTO  800 

ELSE  I F  (DABS(X( 3) )  .GT.  1.012)  THEN 
Cthe  root  has  yet  to  be  identified 
COUNTER =C0UNTER*1 
xc i )=y (2> 

GOTO  503 
ELSE 


C  23.567* 


ROUTINE  TO  FIND  CUB ! C  ROOT 


k4  700  THETA(1)=(3.DO*OSORT(3)*a)/(2.DO*P*DSQRT(  P)) 

'.5  THE  T  A( 1 )=(DAC0S(TMETA(1 ) )  )/3 .00 

1.6  Th£TA(2)=ThETA(1  )*((2.D0*PI >/3 .DO) 

1.7  THETA(3)=THETA(1)  <(2.D0*PI )/3.D0> 

’.8  £=OSORT(( -4.D0*P)/3.D0) 

1.9  K(1)=E*OCOS(TMETA(D) 

180  X(2)=E*DCOS(THETA(2>) 

131  X( 3 )=E*OCOS( THE  TA( 3 ) ) 

132  F(1)=X(1)  <B/(3.D0*C)) 

153  F(2)=X<2)  (B/(3.D0*C)> 

13*.  F(3)=X(3)  (8/ (3.00*0) 

153  00  10  i=1,J 

136  UR  I TE ( 1 ,210)  1,  F(i) 

157  WRITE(*,210)  i,  F(j) 

158  10  CONTINUE 

159  GOTO  900 

160  C234567****** *************************************** ************** 

161  800  WRI TEC  1 ,8101  F(1),  X(3l 

162  URITE(*,810)  F(1),  X ( 3 > 

163  810  FORMAT ( IX, 1  The  root  is  ',F30.15./,‘  The  error  is  • , f 30 . 1 5 > 

164  C23456 7********* ******** . *•••*•*•******»**••»*•*•»***#•*••■ 

165  900  STOP 

166  ENO 
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1  C234567** ******* . *•****««**» . 

2  C  CAPMEGA  is  a  program  that  finds  the  root  to  delta  cap 

3  C  omega,  the  equation  delta  cap  omega  being  defined  in 

4  C  The  Motion  of  a  satellite  in  an  Ax i  •  synmetr  ic  Gravi tat ional 

5  C  Field,  by  R .  H.  Merson  .  As  found  in  the  Geophysical  Journal 

6  C  Vol  4,  1961,  p.17.  This  program  finds  the  roots  of  delta 

7  C  omega  as  a  function  of  f,  where  f=(sin  i)**2.  The  primary 

8  C  equation  is  :  delta  cap  omega  =  Af**2  ♦  Bf  ♦  C  =  0.  Where 

9  C  A,B,and  C  are  constants  depending  on  the  semi  major  axis  and 

10  C  the  Eccentricity 

11  C 

12  C  WRITTEN  BY  J.W.  FOISTER,  III 

13  C  DATE  AUG  21  1987 

14  C 

15  C  J(1)  = 

16  C  J ( 2 )  =  C20  coefficient 

17  C  J(3)  = 

18  C  J(4)  =  C40  coefficient 

19  C  J(5)  = 

20  C  J(6)  =  C60  coefficient 

21  c234567**********«**»********— ************************************ 

22  SNOFLOATCALLS 

23  PROGRAM  CAPMEGA 

24  C 

25  C 

26  C234567**********  DEFINE  THE  VARIABLES  *************************** 

27  REAL*8  J<6) , LAT , A, B,C ,0 , E , F(3) ,EC,DLTA, P,Q, R , X(3> , Y , THETAC3 > , 

28  1  SEM1,P1,W(3),TP,AHLD,RP 

29  C 

30  C 

31  C234567 

32  INTEGER  i.k, COUNTER 

33  C 

34  PI=4.D0*DATAN(1 .00) 

35  R=1.D0 

36  C234567 

37  J<2)=-0. 1960454460-02 

38  J(4)=0. 1889436780-04 

39  J(6)=0. 1340756980-05 

40  C 

41  OPEN  ( 1 ,  F  HE* 1 LPT1 1  ) 

42  200  FORMAT ( IX , ' f  *  ' ,F30.15) 

43  210  FORMAT  ( IX, ' f ' , i 1 , 1  =  \F30.15) 

44  C 

45  C234567*******  INPUT  THE  DATA  ***»»*»**»***»***»****»******»**** 


46 

WRITEC*, 100) 

47 

100 

FORMAT(1X, 'Semi  Major  axis  equals.. 

.  in  KMs _ (F30.15) _ 

48 

READ(*, ' (F30. 15) ' )  SEMI 

49 

WR I TE<*, 1 10) 

50 

11) 

51 

READ( *, 1 ( F30. 15) 1 )  EC 

52  C234567  find  the  Radius  of  the  Peri  apsis 

53  RP=SEMI*(1  EC) 

54  WR!TE(*,105)  SEMl.RP 

55  WR I TE( 1 , 105 )  SEMl.RP 

56  105  FORMAT </, IX, 1  The  Semi  Major  axis  in  Km  is  ',F30.15,/, 

57  I  '  The  Radius  of  Periapsis  in  Km  is  • , F30 .15) 

58  C234567  convert  from  KM  to  Ou's  this  is  for  a  Mars  orbit 

59  SEM 1 =SEM l /3393 .400 

60  C234567 

61  AHID=( 3.00/2.00) 

62  TP=(2.D0*PI )*( SEMI **AHL0 ) 

63  TP=TP*( 15. 9197403800) 

64  C 

65  WR!TE(*,120)  SEMI, TP, EC 

66  WR I TE{ 1 , 1 20)  SEMI, TP, EC 

67  120  F0RMAT(1X, 'For  the  Semi  Major  axis  =  '.F30.15, 

68  1  1  the  period  (in  minutes)  is  ■ , F30 . 15 , / , 

69  2  ■  and  Eccentricity  *  '.F30.15) 


70  C23-567*********  CALCULATE  THE  VALUES  OF  THE  COEFFICIENTS  ********* 

71  LAT  =  SEMI*d  ,00  EC**2) 

72  C 

73  J  ( 1 )=J(4)*( (R/L AT )**4) 

74  J(3)=J(6)*((R/LAT)**6) 

75  J(5)=(J(2)**2)*((R/LAT)**4) 

76  C234567 

77  A  =  ((  51975.00/1024. 00)*EC«*4 )  ( (1 7325 .00/128. DC ;' '>.C**2) 

78  1  (3465.00/128.00) 

79  A=J(3)*A 

80  C 

81  8=J(3)*(((141  75 .00/256 .00) *£C**4  )+({4725 .D0/32.D0)*EC**2}+ 

82  1(945.00/32. 00) )*(J(1)*((( -31 5. 00/32.00 )*EC**2> - (105 .D0/16.u0) 

83  2  ))+(J(5)*(((-45. 00/32. 00)*EC**2) - < 15. 00/2. DO))) 

84  C 

85  C=J(3)*(((- 1575. 00/128. 00)*EC**4)<(-525. DO/16. D0)*EC**2> 

86  1  (105. 00/16. 00))*(J(1)«(((45. 00/8. 00)*EC**2)*( 15 . 00/4. DO))) 

87  2  ♦(J<5)*<((-3. 00/8. 00)*EC**2)*(9. DO/8. 00)>) 

88  C=C*(J(2)*((R/LAT)*«2))*( -3.00/2.00) 

89  C234567 

90  U(1)=(B**2)-(4.D0*A*C) 

91  IF  (W(1)  .LT.  0.0)  THEN 

92  UR1TE(*,300) 

93  WR!TE(1,300) 

94  300  FORMAT ( IX , ' THE  ROOTS  ARE  IMAGINARY!1) 

95  W(2)=-B/(2.D0*A) 

96  U(3)=W(1)/(2.D0*A> 

97  HI RITE(*,310)  U(2),W(3) 

98  UR  I TE ( 1 ,310)  U(2),U(3) 

99  310  F0RMAT(1X, 'ROOT  1  =  ',£30.15,'  +  i  ',F30.15> 

100  WRITE!*, 320)  U(2) ,U(3) 

101  WR I TE( 1 , 320)  U(2),W(3) 

102  320  F0RMAT(1X, 'ROOT  2  =  '.F30.15,'  •  i  ',F30.15) 

103  GOTO  900 

104  ELSE 

105  END  IF 

106  C 

107  Fd)»(-8*DSORT(Ud)))/(2.DO*A) 

108  F ( 2 )=(  8  OSQRT(U(1)))/(2.D0*A) 

109  URITE(*,330)  F(1),F(2) 

110  WRITEd  ,330)  F(  1 ) ,  F(2) 

111  330  FORMAT! IX, 'ROOT  1  =  ',F30.15,/,'  ROOT  2  =  ',F30.15) 

112  C234567 

113  900  STOP 


s 


Program  Dsemi 


C234567 . . 

C  This  program  solves  for  the  change  in  the  semi  major  axis, 

C  and  the  change  in  the  eccentricity  over  one  orbital  period 

C  yith  the  change  due  solely  to  air  drag.  Equation  4.14  and 

C  equation  4.11  from  Desmond  King-Hele's  book,  THEORY  OF 

C  SATELLITE  ORBITS  IN  AN  ATMOSPHERE  are  used  as  the 
C  expression  of  change  in  semi  major  axis,  and  eccentricity 

C  These  equations  include  integrals,  and  therefore  requrires 

C  an  integration.  A  standard  8th  order  Gaussian -Legendre 
C  quadrature  method  is  used  to  perform  this  integration. 

C  The  interval  of  integration  (from  0  to  2  pi )  wi i l  be  broken 
C  up  into  four  intervals  (from  0  to  pi/2,  from  pi/2  to  pi 

C  from  pi  to  Spi/2,  from  3pi/2  to  2pi )  and  the  Gaussian  - 

C  Legendre  quadrature  will  be  applied  to  each  interval.  This 
C  will  improve  the  accuracy  of  this  routine. 

C 

CWritten  by  J.  W.  Foister,  III 
C 


19  CDate 


22  Sept.  1987 


C 

C234567************ ****** ************************************ 

SNOFLOATCALLS 

PROGRAM  DSEMI 
C 

REAL*8  A,AK(4),EK(4),PI , RHO, RHOO .DELTA, F , INCL,RPO,VPO, 

1  S,CD,M,EC,H,MU,U,DA,DE,E,YE,XE, I  NT 


27  C234567*** 


C  EC 

C  H 

C  MU 

C 

c  u 

C  DA 

C  OE 

C  YE 

C  XE 

C  DE 

C  INT 

C 

C234567** 

C 


•DEFINITION  OF  VARIABLES****************************** 

=  intial  value  of  the  semi  major  axis  (km) 

=  the  Gaussian  quadrature  coefficients 
=  the  values  of  E,  the  independent  varialbe  where 
f(E)  =  0. 

*  the  irrational  nurber  pie. 

=  the  density  of  the  atmosphere  in  (kg)/(km**3) 

=  the  reference  density  of  the  atmosphere 

*  the  surface  area  of  the  satellite  times  the  coeff. 
of  drag  divided  by  the  mass  of  the  satellite,  with 
the  entire  quanity  multiplied  by  a  correction 
factor  that  converts  the  velocity  of  the  satellite 
wrt  the  center  of  the  planet,  to  a  velocity  wrt  the 
atmosphere  in  which  it  is  moving.  (F) 

=  the  reference  periapsis  (km) 

=  the  conversion  factor  that  changes  the  velocity 
term  from  one  relative  to  the  center  of  the  planet 
to  one  relative  to  the  atmosphere 
=  the  angle  of  inclination  (radians) 

=  the  velocity  of  the  satellite  at  RPO  (km)/(sec) 

=  the  reference  area  of  the  satellite  (km**2) 

=  the  coefficient  of  drag 
i  the  mass  of  the  satellite  (kg) 

=  the  eccentricity 
=  the  scale  height  (km) 

=  the  gravitational  constant  times  the  radius  of  the 
planet  (km»*3)/(sec**2) 

=  the  rotational  velocity  of  the  planet  (rad/sec) 

=  the  change  in  A 
=  the  change  in  Eccentricity 
=  the  integrand  for  change  in  A  equation 
=  the  integrand  for  change  in  Eccentricity  equation 
=  the  change  in  eccentricity  over  one  orbital  period 
=  the  constant  need  inorder  to  perform  the  needed 
change  in  varialbe  required  by  the  interval 


INTEGER  i,j,k 


DATA  Ek/. 9602898600, . 7966664800, . 5255324 IDO, . 1834346400/ 
DATA  AK/. 1012285400, .2223810300, .3137066500,-3626837800/ 
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71  OPENd.FILE^LPTI') 

72  C 

73  C234567*************DEFINE  THE  CON ST AMTS* ******** 

74  PI=4.D0*0ATAN(1 .DO) 

75  RHO0=3.30-3 

76  S=1 -D-  5 

77  CO=2.DO 

78  M=1000.00 

79  H= 14. 1386704900 

80  RPO=3593.4DO 

81  MU=4.2828287D4 

82  VPO=DSQRT(MU/RPO) 

83  W=<4. 0612498030 • 3>*(PI/180. DO) 

84  !NCL=(45.D0)*(PI/180.D0) 

85  F  =  (1 .DO-<RPO/VPO)*W*DCOS(INCL))**2 

86  A=5133.428571D0 

87  EC=.300 

88  DELTAIC F*S*CO )/M 

89  WR I TE(* , 500)  PI  ,RHOO,S,CO,M,H,RPO,MU 

90  UR! TE( 1 ,500)  PI , RHOO, S, CD ,M, H, RPO.MU 

91  500  FORMAT ( 1 X , 1  PI  =  * , P30 .15,'  RHO0=  • , F30.15,/ 

92  1  S=  ' , F30. 15, 1  CD=  '.F30.15,/, 

93  2  >  M=  >,£30.15/  H=  1 , F30. 15,/, 

94  3  '  RPO=',F30.15, >  MU=  >,F30.15) 

95  WRI TEC*, 600)  VPO,W, I NCL , F , A, EC.DELTA 

96  WR I TE( 1 ,600)  VPO.W, INCL , F , A, EC .DELTA 

97  600  FORMAT < 1 X, >VPO=  • , F30 .15,'  W=  • , F30 . 15 , / , 

98  1  1  I NCL— ‘ .F30.15, '  F=  ' ,F30. 15,/, 

99  2  •  A  =  1 , F30. 15, 1  EC=> ,F30.15,/, 

100  3  ■  DELTA  =  1 ,F30. 15) 

101  C234567****************************************** 

102  DA=O.DO 

103  DE=0 .00 

104  C234567*********determine  the  interval*********** 

105  00  10  1=1,4 

106  IF  (i  .EQ.  1)  THEN 

107  INT=1.D0 

108  ELSE  I F  (i  .EQ.  2)  THEN 

109  INT=3.00 

110  ELSE  I F  (i  .EQ.  3)  THEN 

111  INT=5 .00 

112  ELSE  I F  <i  .EQ.  4)  THEN 

113  INT=7.D0 

114  ELSE 

115  ENDIF 


116 

C234567********start  the  quadrature******************* 

117 

00  20  J=1 ,4 

118 

DO  30  k=1 ,2 

119 

IF  <k  .EQ.  1)  THEN 

120 

E=(PI/4.00)*<EK(j)*INT) 

121 

ELSE I F  (lc  .EQ.  2)  THEN 

122 

E=(P1/4.D0)*<  EK< j)*INT) 

123 

ELSE 

124 

ENDIF 

125 

RHO=RHOO*DEXP<(< -A*EC)/H)*(1 .DO-DCOS(E))) 

126 

YE=(1.00*EC*DCOS(E))**(1.5) 

127 

YE=(YE/OSORT<1 .DO  - EC*OCOS( E ) ) >*RHO 

128 

XE=RHO*( 1 .DO- EC**2)*DC0S(E ) 

129 

XE=XE»DSQRT((1.00+EC*0COS(E))/(1.D0-EC*DCOS(E 

130 

DA=( AX( j )*YE )*0A 

131 
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